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SUMMARY 


A time-domain finite element method is developed for optimal control prob- 
lems. The theory derived is general enough to handle a large class of problems 
including optimal control problems that are continuous in the states and controls, 
problems with discontinuities in the states and/or system equations, problems with 
control inequality constraints, problems with state inequality constraints, or prob- 
lems involving any combination of the above. The theory is developed in such a way 
that no numerical quadrature is necessary regardless of the degree of nonlinearity 
in the equations. Also, the same shape functions may be employed for every prob- 
lem because all strong boundary conditions are transformed into natural or weak 
boundary conditions. In addition, the resulting nonlinear algebraic equations are 
very sparse. Use of sparse matrix solvers allows for the rapid and accurate solution 
of very difficult optimization problems. The formulation is applied to launch- vehicle 
trajectory optimization problems, and results show that real-time optimal guidance 
is realizable with this method. Finally, a general problem solving environment is 
created for solving a large class of optimal control problems. The algorithm uses 
both FORTRAN and a symbolic computation program to solve problems with a 
minimum of user interaction. The use of symbolic computation eliminates the need 
for user-written subroutines which greatly reduces the setup time for solving prob- 
lems. 
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CHAPTER 1 


INTRODUCTION 


The calculus of variations was born at the very end of the seventeenth century 
through the work of such great mathematicians as Newton and Leibniz. With 
the birth of the calculus of variations came optimal control theory for continuous 
systems. The newly discovered methods of differential calculus were used at once to 
solve many important and practical maximum and minimum problems. The first 
optimization problem solved by the calculus of variations was set up and solved 
by Newton in 1686. The problem, of interest to aerospace engineers even today, 
was to choose a nose shape for minimum drag in hypersonic flow [1]. Another 
optimization problem was presented by John Bernoulli in 1696. Bernoulli posed 
the classical brachistochrone problem defined as: Among all lines connecting two 
given points, find the curve traversed in the shortest time by a material body under 
the influence of gravity [2]. This problem was solved independently by John and 
James Bernoulli, L’Hopital, Leibniz, and Newton. For years, work was done on 
solving separate variational problems. However, the tremendous intellectual feat 
of creating a single method for solving variational problems belongs to the great 
Swiss mathematician Leonhard Euler. At the astonishingly young age of 25, Euler 
published his work “General Solution of the Isoperimetric Problem Taken in Its 
Most General Sense” [3]. 

Since the birth of the calculus of variations, optimization problems have been 
a topic of research. The optimal control problem of interest in this thesis may 
be described as follows. Consider a system that is completely defined by a finite 



number of states, i.e., quantities that describe the current status of the system. 
The status of the system is determined by a set of first-order ordinary differential 
equations. The states are influenced by a finite number of control variables. The 
optimization problem is to choose the control variables to satisfy the given boundary 
conditions while minimizing (or maximizing) a given performance index, or cost 
functional. Use of the calculus of variations results in a multi-point boundary- value 
problem. Unfortunately, not many analytical solutions to these types of optimal 
control problems have been found beyond that which the great minds of the 17 th 
and 18 tft centuries solved. However, the appearance of practical, high-speed digital 
computers in the 1950’s revolutionized the field of optimal control. Computers and 
numerous numerical methods are now the tools for dealing with the nonlinear and 
complex systems of today. Still, with the ever-increasing attention given both space 
exploration and space travel, even more reliable and efficient numerical methods 
are required. This thesis deals with a new type of numerical method based on finite 
elements in time for solving optimal control problems. Particular emphasis is given 
to the computation of optimal trajectories for advanced launch vehicles. 


Background 


If the United States is to maintain its position as a world leader among space- 
faring nations, then cheaper and more reliable means for transporting people and 
cargo to and from space must be developed. The current Space Shuttle, as tech- 
nically successful as it is, will not meet all the future needs of the United States. 
Studies indicate that the projected transportation needs will be best served by a 
mix of expendable and reusable vehicles. Specifically, the functions of one-way cargo 
transport to orbit and two-way passenger transport should be separated [4]. 

The Air Force turned to industry in May of 1987 to help meet the goals of the 
advanced launch system (ALS) program. The goals of the ALS program are to (1) 
place large payloads (in excess of 100,000 pounds) into low Earth orbit and at an 
order of magnitude lower cost per pound, and (2) make space launch operations, 
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including the manufacturing, processing, and actual launch of booster vehicles, 
significantly more routine compared to present methods and procedures [5]. 

Currently, an extensive amount of ground support (typically weeks) is required 
to prepare the guidance system of the Space Shuttle for launch. To meet the oper- 
ational requirements of the ALS program, ground support for pre-mission activities 
must be drastically reduced. On-board algorithms must maximize system perfor- 
mance as measured by autonomy, mission flexibility, in-flight adaptability, reliabil- 
ity, accuracy, and payload capacity. For real-time trajectory optimization to be 
realizable, the algorithms must be computationally efficient, robust, self-starting, 
and capable of functioning independently of ground control. Furthermore, the algo- 
rithms must be designed with the anticipation that the launch vehicle will undergo 
evolutionary growth [6]. 

On-board, real-time trajectory optimization algorithms are required to meet 
the needs of the ALS program. Such algorithms promise to (1) reduce the cost 
of designing flight profiles, (2) reduce the time required to respond to changed 
payload or mission requirements, and (3) improve vehicle performance. The diverse 
mission requirements of a general-purpose launch vehicle require new approaches to 
trajectory optimization. This work, conducted under NASA Grant NAG-1-939 of 
which Dr. Daniel D. Moerder is the technical monitor, was chiefly concerned with 
the development of a method for calculating optimal trajectories of these advanced 
launch vehicles. 


1.2. Previous Work 

There are scores of methods available for solving optimal control problems. 
The method chosen is dependent upon, among other things, the type of problem to 
be solved and the resources (in terms of both software and money) available to the 
user. This section will give a brief outline of some of the methods now being used 
to solve optimal control problems. 
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1.2.1. The LQ Problem 


One type of problem considered frequently is one in which the system is defined 
by a set of linear differential equations while the performance index is a quadratic 
functional. This is the so-called LQ (linear-quadratic) problem. Over the past 
twenty years, a new technique has been established for the solution of LQ prob- 
lems using orthogonal functions. The main characteristic of this technique is that 
the differential equations involved in the problem are reduced to algebraic equa- 
tions, thus greatly simplifying the problem solution. The technique calls for the 
differential equations to first be converted to integral equations. Subsequently, the 
various unknowns involved in the integral equation are approximated by truncated 
orthogonal series. The key idea of this technique is to derive an integral operational 
matrix to convert integral expressions into algebraic equations. The form of the 
operational matrix is dependent upon the choice of orthogonal functions used. Var- 
ious functions have been used to parameterize the system including Walsh functions 
[7, 8], block-pulse functions [9, 10], Chebychev functions [11], Hermite series [12], 
polynomial series basis vectors [13, 14], and Legendre polynomials [15]. 

Another approach to LQ problems is presented in [16]. The method transforms 
the canonical equations to a set of algebraic equations and allows approximating 
functions that need not satisfy the initial conditions a priori. This enlarges the 
space from which the approximating functions can be chosen. Furthermore, a La- 
grange multiplier technique is used to enforce the terminal conditions on the states. 
Orthogonal polynomials are then used to solve the LQ problem. The idea of setting 
up approximating functions that do not need to satisfy boundary conditions is one 
of the key ideas in this current work, as will be seen in Chapter 2. 

1.2.2. The Nonlinear Problem 

Methods available for the solution of optimal control problems generally fall 
into two distinct categories: direct and indirect. Direct techniques seek to directly 
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minimize the performance index, or cost functional, by prudent choices of the free 
parameters in the system. Indirect techniques, on the other hand, seek to minimize 
the performance index indirectly by satisfying the first-order necessary conditions 
for optimality as established from the calculus of variations. 

The direct approach to the solution of optimal control problems first requires 
parameterization of the control and state time histories. The choice of parameteri- 
zation schemes is not unique and success of the direct methods has been achieved 
using Hermite polynomials [17], Chebychev polynomials [18, 19], single-term Walsh 
series [20], splines [21], and the like. 

Once the parameterization scheme is chosen, a parameter optimization algo- 
rithm is then used to improve the initial guess of the free parameters. These algo- 
rithms are in common use today and include variable metric techniques or quasi- 
newton methods [22] and variations on gradient methods. Gradient methods [1, 
23] were developed to surmount the “initial guess” difficulty associated with other 
methods such as Newton-Raphson. They are characterized by iterative algorithms 
for improving estimates of the control histories, in order to come closer to satis- 
fying the optimality conditions and the boundary conditions. First-order gradient 
methods usually show rapid improvements when sufficiently far from the optimal 
solution. However, the rate of convergence drastically decreases in the neighbor- 
hood of the solution. Second-order gradient methods have excellent convergence 
characteristics near the optimal solution, similar to a Newton-Raphson method. 
Conjugate gradient methods are very powerful because they combine the first-order 
and second-order gradient methods. Ref. [24] contains a thorough description of 
the gradient method and many other algorithmic methods in optimal control. It is 
noted that direct methods have been successfully used to solve trajectory optimiza- 
tion problems [17, 18, 25, 26]. 

The indirect approach to the solution of optimal control problems attempts to 
satisfy the necessary conditions of optimality as derived from the calculus of varia- 
tions. These conditions result in multi-point boundary-value problems. Analytical 
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solutions to such problems are generally unobtainable except for the simplest of 
problems. Therefore, numerical methods are usually employed. 

The two main techniques for solving nonlinear multi-point boundary-value 
problems are shooting methods and quasilinearization methods. Shooting meth- 
ods [27, 28, 29] are frequently used and can be described as follows: The initial 
conditions and the differential equations are satisfied at each stage of the process 
while the final conditions are sacrificed somewhat. A nominal solution is generated 
by guessing the missing initial conditions and forward integrating the differential 
equations. The intent is to reduce the error in the final conditions at each itera- 
tion. Quasilinearization techniques [1, 30] involve choosing nominal functions for 
the states and costates that satisfy as many of the boundary conditions as possible. 
The control is then found by using the optimality conditions. The system equations 
and costate equations are then linearized about the nominal values and a succession 
of nonhomogeneous, linear two-point boundary value problems are solved to modify 
the solution until the desired accuracy is obtained. Other indirect techniques in- 
clude the method of adjoints, Newton- Raphson methods, and continuation methods 
[31, 32, 33] 

1.2.3. Finite Element Methods 

Finite element methods, which include the Rayleigh- Ritz and Galerkin methods 
[21], as well as the method of collocation have been used to solve optimal control 
problems [34]. Ref. [35] appears to be one of the first papers using finite elements 
to solve optimal control problems. Therein, the authors considered the application 
of a modified Ritz-Trefftz direct method to the so-called state regulator control 
problem. This is an LQ problem and the method leads to an approximation of the 
performance index of order h 7 where h is the time step involved. The modified Ritz- 
Trefftz method used in [35] was later extended to include problems with terminal 
conditions on the states [36]. Other examples of the use of the Ritz method can be 
found in [37] and [38]. An application of the collocation method is found in [17]. 
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Only a few of these papers are listed here because these methods, although very 
accurate and useful, suffer from the same computational problems. These methods 
require that the approximating functions satisfy all the strong boundary conditions. 
Thus, in practice, certain equations are eliminated depending on the boundary 
conditions present for a particular problem. Another drawback of these methods is 
that numerical quadrature is required. This can introduce error and greatly increase 
the computational effort required to solve a problem. The finite element method 
described in this thesis avoids these two pitfalls, yielding a computationally efficient 
and versatile algorithm. 

1.2.4. General Programs 

In closing this section, a few of the commercially available programs for solving 
optimal control problems are mentioned. The first three programs mentioned are 
general-purpose problem solvers, whereas the last two are particularly designed to 
optimize point-mass trajectories. 

The Chebychev Trajectory Optimization Program (CTOP) has been found to 
be useful in a wide variety of practical applications [18]. This program uses a direct 
technique for solving problems and parameterizes the functions using Chebychev 
polynomials. Penalty functions are used to enforce the equations of motion and 
path constraints. The Nonlinear Programming for Direct Optimization of Trajecto- 
ries (NPDOT) uses piecewise polynomials and collocation to satisfy the differential 
equations. Results presented in [17] show that NPDOT runs much more quickly 
than does CTOP. A FORTRAN program called MISER (the origin of the acronym 
is unknown) is presented in [39]. This general-purpose software utilizes a unified 
computational approach to solve a wide range of optimal control problems subject 
to general constraints. This program appears very useful; however, a substantial 
amount of user programming is involved. Specifically, the user must provide a 
series of FORTRAN subroutines which evaluate the right-hand side of the state 
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and costate equations. Furthermore, the user must transform the problem into the 
canonical form described by the authors. 

POST, or Program To Optimize Simulated Trajectories, provides the capability 
to target and optimize point mass trajectories for a powered or unpowered vehicle 
operating near a rotating oblate planet [40]. POST offers the solution to a wide 
range of flight problems including aircraft performance, orbital maneuvers, and 
injection into orbit. The user can select the optimization variable, the dependent 
variables, and the independent variables from a list of more than 400 program 
variables, POST is also operational on several computer systems. Another much- 
used program is OTIS, Optimal Trajectories by Implicit Simulation. OTIS is a 
three degree of freedom (point mass) simulation program for multiple vehicles [17]. 
The user can simulate a wide variety of vehicles such as aircraft, missiles, re-entry 
vehicles, and hypervelocity vehicles. The methods used were chosen to improve 
speed, convergence and applicability of OTIS over existing performance programs. 
Both POST and OTIS are very reliable and accurate programs, but are limited in 
scope as compared to the three programs listed above. 

1.3. Present Approach 

This thesis describes in detail a time-domain finite element approach for solving 
optimal control problems. The so-called weak principle for optimal control problems 
is based on Hamilton’s principle, which has traditionally been used in analytical me- 
chanics as a method of obtaining the equations of motion for dynamical systems. 
Bailey [41] followed by several others [42, 43, 44] obtained direct solutions to dy- 
namics problems using a form of Hamilton’s principle known as the law of varying 
action, thus opening the door for its use in computational mechanics. 

More recently, it has been shown that expression of Hamilton’s law as a weak 
form (commonly referred to as Hamilton’s weak principle or HWP) provides a pow- 
erful alternative to numerical solution of ordinary differential equations in the time 
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domain [45, 46]. The accuracy of the time-marching procedure derived in [45, 46] 
is competitive with standard ordinary differential equation solvers. Further com- 
putational advantages were obtained in so-called mixed formulations of HWP in 
which the generalized coordinates and momenta appear as independent unknowns 
[47]. Therein, an unconditionally stable algorithm emerges for the linear oscilla- 
tor with exact element quadrature. HWP also has shown to be an ideal tool for 
obtaining periodic solutions for autonomous systems, as well as finding the corre- 
sponding transition matrix for perturbations about the periodic solution [48]. These 
are complex two-point boundary value problems; its utility for these problems and 
its superior performance in mixed form strongly suggest that it could be used in 
optimal control problems. 

Chapter 2 develops the weak principle for optimal control theory for problems 
in which the states and controls are continuous. In Chapter 3, the theory developed 
in Chapter 2 is tested on a single-stage rocket trajectory optimization problem. The 
weak principle is then extended in Chapter 4 to handle problems with discontinuities 
in the states and system equations. A realistic two-stage rocket problem is then 
solved in Chapter 5. The weak principle is further extended in Chapters 6 and 7 
with the inclusion of control and state inequality constraints. A final demonstration 
of the trajectory optimization capabilities of the weak principle is demonstrated in 
Chapter 8. 

Chapter 9 of this thesis gives a brief study of error estimates for the weak 
principle. The subject of generating initial guesses to solve the discretized equations 
is dealt with in Chapter 10. Chapter 11 is the main contribution of the work. It 
describes a general code for the solution of optimal control problems. Conclusions 
and future research are discussed in Chapter 12. 

There are three appendices to the thesis. Appendix A discusses the solution 
of dynamics problems using Hamilton’s weak principle. Appendix B deals with the 
solution of initial- value ordinary differential equations by using the weak principle. 
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Finally, Appendix C describes how simple beam problems can be cast in the form 
of optimal control problems and solved using the weak principle. 
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CHAPTER 2 


WEAK PRINCIPLE FOR OPTIMAL CONTROL 


It is desired to develop a solution strategy for optimal control problems based on 
finite elements in time. Finite elements have been used in the past to solve optimal 
control problems and two-point boundary-value problems in general (see [21] and 
[35-38]); however, these methods all require numerical quadrature. In addition, a 
different choice of shape functions must be made for each problem depending on 
the strong boundary conditions to be enforced. These two obstacles are overcome 
with the weak principle derived below. 

A weak principle based on the variation of the performance index will be for- 
mulated [49, 50]. When deriving this formulation, two things must be remem- 
bered. First, the resulting formulation must satisfy the Euler-Lagrange equations 
and boundary conditions that have already been established in optimal control the- 
ory [1]. Second, in an attempt to make the solution scheme as general as possible 
all strong boundary conditions will be transformed into natural or weak boundary 
conditions. 

The boundary conditions axe all cast in the form of weak boundary conditions 
so that the shape functions used for the test functions can be chosen from a less 
restrictive class of functions. For example, if there is a strong boundary condition 
on one of the states at the initial time (*.e., an initial condition) then the shape 
function chosen for the variation of that state must equal zero at the initial time 
[51]. It would be advantageous if one could choose the same shape functions for 
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every optimal control problem. This is possible if there are no strong boundary 
conditions that must be satisfied by the shape functions. 

The idea of transforming strong boundary conditions to natural boundary con- 
ditions [52] revolves around adjoining a constraint equation to the performance 
index with an unknown Lagrange multiplier. The variation of the performance in- 
dex is then taken in a straightforward manner. Through appropriate integration 
by parts, one may show that the Euler-Lagrange equations are identical to those 
derived in classical textbooks [1] and that the boundary conditions are the same, 
only stated weakly instead of strongly. 

As is shown below, the weak principle for optimal control reduces the necessary 
conditions for optimality to a set of nonlinear algebraic equations. These algebraic 
equations can be derived prior to specifying the problem to be solved. It is this 
feature in particular that makes the weak principle so powerful. 

2.1. General Development 

Consider a system defined by a set of n states x and a set of m controls u. 
Furthermore, let the system be governed by a set of state equations of the form 
x = f(x,u,t). In this chapter, the class of problems is restricted to those where x, 
u , and / are continuous. We may denote elements of the performance index, J 0 , 
with an integrand L(x,u,t ) and discrete functions of the states and time (f>\x(t),t] 
defined only at the initial and final times t o and tf. In addition, any constraints 
imposed on the states and time at the initial and final times may be placed in sets 
of functions ^>[x(t), <]. These constraints may be adjoined to the performance index 
by discrete Lagrange multipliers v defined at t 0 and tf. Finally, we may adjoin the 
state equations to the performance index with a set of Lagrange multiplier functions 
A(t) which will be referred to as costates. For variable tf, this yields a performance 
index of the form 
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where $ = cf)[x(t),t] + u T ip[x(t),t]. The constraints to be adjoined to Jo above are 
simply that the states be continuous at the initial and final times. Introducing 


and 


®|i 0 = lim x{t) and :r| t/ = lim x(t) (2-1-2) 

t-tf 7 


x 0 = x\t 0 = x(to) and x f = x\ tf = x(t f ) (2.1-3) 

continuity is weakly enforced by adjoining a T (x — i)|(„ to Jo where a is a set of dis- 
crete unknown Lagrange multipliers defined only at to and tf. The new performance 
index is 


J = f [£(*.«,() + A r (/ - ±)] dt + $11; + c T (x - x)|Jj (2.1-4) 

Jt 0 

To derive the weak principle, it is necessary to determine dJ , the first variation 
of J. Denoting with 5x(tf) and 8\ (tf) the variations of x and A at t = tf when 
holding t f fixed, and letting dx(tf) and d\(t /) be the variations of x and A at t = tf 
when tf is allowed to be free, then the variations at t = tf can be expressed by the 
linear equations (see [1] or [2]) 


8x(tf) — dx(tf) - x\ t j dtf and 8\(tf) = d\(tf) - A| t/ dtf (2.1-5) 


The first variation of J is 




( 2 . 1 - 6 ) 


A necessary condition for an extremal of J is that the first variation be zero. 
Also, the admissible variations of the states must be continuous at the initial and 
final times and therefore ( dx — = 0. For notational convenience we will define 



dx 


to 


and 




(2.1-7) 


As an aside, to ensure that no necessary conditions have been altered, the 8x term 
will be integrated by parts. This results in 



+ ixT K[ + 6 ° T(x - ^i;:= ° 


( 2 . 1 - 8 ) 


Using Eq. (2.1-5) and noting that <5ar| to = dx | <0 since is fixed, then the above 
equation can be simplified to the following. 
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(2.1-9) 



The Euler- Lagrange equations from the above will now be compared with the well- 
known optimal control equations presented in [1]. The coefficients of 8 , 8x T , and 
8u t in the integrand, when set equal to zero, correspond to Eqs. (2.8.15 - 2.8.17) 
from [1]. There are also four trailing terms in Eq. (2.1-9) from which the boundary 
conditions of [1] can be determined. Namely, the requirement for the coefficient of 
dtf to vanish is equivalent to Eq. (2.8.20). The requirement for the coefficient of 8v T 
to vanish at t = </ yields Eq. (2.8.21). The requirement for the coefficient of dx T 
to vanish at t = if shows that the value of A|^ equals A|( 7 as given in Eq. (2.1-7), 
which corresponds to Eq. (2.8.19). Finally, the requirement for the coefficient of 
8<x t to vanish at t = to requires the value of x\ to to equal x|t 0 , in accordance with 
Eq. (2.8.18). 

Three additional boundary conditions are present in the above formulation. 
One is the requirement for the coefficient of 8a T to vanish at t = tf which demands 
that the value of x\ t/ equal x\ t/ . The second is the requirement for the coefficient 
of dx T to vanish at t = to which demands that the value of A|< 0 equal A|t 0 as given 
in Eq. (2.1-7). These two conditions enforce continuity of the states at the final 
time and continuity of the costates at the initial time. The third and last boundary 
condition is the requirement for the coefficient of 8u T at t = to to vanish which 
demands that ip[x(t 0 ), t 0 ] = 0. Again, all boundary conditions were cast in the 
form of natural boundary conditions so that the shape functions chosen for 8x and 
8A will not have to satisfy any particular boundary conditions. 
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Having satisfied the requirement that none of the fundamental equations are 
altered, the weak formulation is now derived from Eq. (2.1-9). By noting that a 
is a Lagrange multiplier whose only restriction is that 8a be independent of dx , 
8u y and dtf , and that a has the units of the costates, we then choose 8a = dX. 
(Note that there is no unique choice for 8a , but this one will lead to a successful 
solution strategy.) Also, the dx and d\ terms can be eliminated from Eq. (2.1-9) 
using Eq. (2.1-5) resulting in 




+ 8u t 


4* dtf 



</ 

*0 


+ &r T (A - A)|^+<5A r (a: - x)|j' 

+ (A — X) T x\ t ^dtf — (x — x) T \\ t ?dtf = 0 


( 2 . 1 - 10 ) 


Finally, the x and A terms are integrated by parts so that no time derivatives 
of i or A appear in the weak formulation for optimal control. This allows for the 
simplest possible shape functions to be chosen which in turn eliminates the need 
for numerical quadrature. The resulting equation is 


f{ w 


8\ T x — 8 x t A + 8x t 
4- 8\ T f + 8u t 

+ dtf ( L 4- A T / 4” 


dL Y + m T 


dx) + (j&J A 


dL\ 
du ) 

dQY 

dt j 


+ 


m 


dt 


( 2 - 1 - 11 ) 


+ 6v T tp\'' 


to 


+ 8x t A ^ — 8\ t x ^ 4- (A — \) T x f dtf — (x — x) T A| t dtf = 0 
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After noting that the last two terms are equal to zero in accordance with the natural 
boundary conditions (see the coefficients of dx and Sot in Eq. 2.1—9), one can discard 
those terms without changing the necessary conditions. Also, we note that for most 
problems, the initial conditions are given for all n states and thus, in accordance 
with Eq. (2.1-7), all the initial costates axe unknown. Therefore, instead of treating 
elements of v at t = to as unknowns and replacing A| to with these unknowns, we 
will instead treat A| to as unknowns and eliminate the Sv\t 0 equations from the weak 
principle. We hasten to point out that the elements of x| to are the initial conditions. 
The final form of the weak principle is then given as 


Jt |^A r x — Sx T X + Sx T 



+ SX T f + Su T 



+ dtf (z + A T / + |£) 



+ Su^\ t +Sx T X\'-SX T x \/ Q - 0 


(2.1-12) 


This is the governing equation for the weak Hamiltonian method for optimal control 
problems of the form specified. It will serve as the basis for the finite element 
discretization described below for constructing of candidate solutions (i.e., solutions 
which satisfy all the necessary conditions). It should be noted that normally one 
will encounter various types of discontinuities in the states and state equations, as 
well as inequality constraints on the controls and states in problems that deal with 
optimal control. These aspects will be treated in Chapters 4, 6 and 7 respectively. 
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2.2. Finite Element Discretization 


Let us break the time interval from to to tf into N elements. The nodal values 
of these elements are i,- for i = 1, . . . ,N + 1 where to = ti and tf = Now, 

define a nondimensional elemental time r as 


t-tj _ t-tj 
t»-f- 1 ti Ati 


( 2 . 2 - 1 ) 


Note in Eq. (2.1-12) that time derivatives of 8x and 8\ are present. However, 
no time derivatives of x and A exist. Therefore, it is possible to implement linear 
shape functions for 8x and 8 A and constant shape functions for x and A within each 
element. The linear shape functions for the virtual states and costates are 


8x = £zj(l — r) + 8xi+\T 
8\ = £A,-(1 — r) + ^A.q-ir 


( 2 . 2 - 2 ) 


For the states and costates, piecewise constant shape functions are taken to be 


smd 



f Xi 

II 

o 


x = \ 

1 - 

if 0 < r < 1 

(2.2-3) 

1 

l Xi+l 

if r = 1 



[A,- 

II 

O 


A = < 

A, 

if 0 < t < 1 

(2.2-4) 


lA i+a 

if r = 1 



It is important to understand that the equalities xi = x(<o),Ai = A(f 0 ), £n+i = 
x(tf), and Aat+i = A (tf) are enforced as natural (weak) boundary conditions. In 
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other words, the hatted values of x and A at the beginning and end of the time 
interval are the discrete values of x and A that are needed in the weak formulation 
of Eq. (2.1-12). This is clarified below in Figure 2.1 where the time line is broken 
into elements and the nodes are labelled appropriately. 



Fig. 2.1: Time line broken into elements and labelling of nodes 


Finally, note that the time derivatives of u and Su do not appear in the formulation. 
Thus, constant shape functions are chosen for both u and the variation of u. These 
shape functions are 


u = Ui 

Su = 6Ui 


(2.2-5) 


Plugging in the shape functions described for x, A, and u , substituting t = 
tj + rAtj, and carrying out the element quadrature over r from 0 to 1 results in 
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N 


Y,\ Sx * 


i= 1 


- An fdf\ T - At, /3t\ r l T f At, 
- AU fdf\ T - At, fdL\ T ] T (_ 

K 2 V»Ji A ' 2 [dxj, w,+l v 




+ SuJ < A ti 


dL 

du 




+ dt } (L + iJf+^ + v Td ' P ' 


dt 


+ Su T ^\ t 


h 


— 8xJ\q + 8 A^x 0 + 8xJj +1 Xf — 8\Jj +1 xf = 0 


( 2 . 2 - 6 ) 


where /, -- f(x = Xi,u = u^), = X(a; = Xi,u = u,). This is the general algebraic 

form of the weak formulation for all optimal control problems of the form specified. 
Note that if the time t does not appear explicitly in the problem formulation then 
all integration is exact and can be done by inspection. If t does appear explicitly, 
then t may be approximated by a constant value over each element (as are x, A, 
and u ) and the integration may still be done by inspection. Note that the elements 
must be assembled over the entire time interval for this two-point boundary-value 
problem. Only the nodal values (the hatted quantities) of the states and costates 
at the initial and final times appear in the algebraic equations. These remaining 
hatted quantities are the discrete values of the states and costates which appear in 
Eq. (2.1-12). 

Eq. (2.2-6) is a set of nonlinear algebraic equations whose size depends on the 
number of elements N . In fact, if tp\ tj is a g x 1 column matrix and there are N 
elements, then there are 2 n(N + 1) (for 8xi and £A;) + mN (for 8u{) + q (for 8u) 
+ 1 (for dtf ) equations and 2n(N + 2) (for Xi, £o> if, -^o, and A/) + mN (for 
u,) + q (for u) + 1 (for t f) unknowns. Therefore, 2n of the An endpoint values 
for the states and costates (xo,Ao ,£/, and A/) must be specified. In general, xo 
(the initial conditions) is known in accordance with physical constraints. Also, A / 
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can be specified in terms of other unknowns with the use of Eq. (2.1-7). Now we 
have the same number of equations as unknowns. These equations may be used for 
any optimal control problem of the form specified. One simply needs to substitute 
the appropriate /, L , <j>, and boundary conditions into Eq. (2.1-12) for a given 
problem. Note that there is never a need to eliminate any of the equations (except 
the dtf equation for fixed-time problems) as is usual in standard finite element 
practice where strong boundary conditions are enforced by the virtual quantities. 

Normally, Eq. (2.2-6) can be solved by expressing the Jacobian explicitly and 
using a Newton-Raphson solution procedure. For the example problems which 
follow, the iteration procedure will converge quickly for a small number of elements 
with a trivial initial guess. Then, the answers obtained for a small number of 
elements can be used to generate initial guesses for a higher number of elements. 
Thus, a large number of elements can be solved with a very efficient run-time on 
the computer. 

Although the nodal values x, and A; for 2 < i < N (on the interior of the time 
interval) do not appear in the algebraic equations, their values can be easily recov- 
ered after the solution is found. This is most easily seen by looking at the following 
ordinary differential equation multiplied by a test (or weighting) function 8X and 
integrated over some time interval where the integral makes sense. A constraint 
to transform the strong boundary conditions to weak ones has been adjoined via a 
discrete multiplier which has been identified as 8\. 

r 8\[f(x,t) -i]dt + 8\(x - x)|j’ = 0 (2.2-7) 

Jti 


After an integration by parts, using the linear shape function for 8 A defined in 
Eq. (2.2-2), using the piecewise constant shape function for x defined in Eq. (2.2- 
3), and substituting r for t as given in Eq. (2.2—1), then the following equation is 
obtained from Eq. (2.2-7). 
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( 2 . 2 - 8 ) 


^Ai(— Xi + + * 1 ) + 8\ 2 (xi + — / - x 2 ) = 0 

With arbitrary £Ai and <5 A 2 , the coefficients must vanish, forming two equations of 
the form 


-X\ + — / + X1 = 0 

At , „ 

xi + —J - x 2 =0 


(2.2-9) 


Now, by subtracting the second equation from the first, it is seen that 


X! 


Xl + X 2 


( 2 . 2 - 10 ) 


or, in words, that the interior value (the bar value) is simply the average of the sur- 
rounding nodal (or hatted) quantities. Once the solution is found, all the midpoint 
values and the end nodal values are known, and thus all other nodal values can be 
recovered by repeatedly using Eq. (2.2-10). In fact, the nodal values are really the 
best approximation to the solution and thus are the only ones plotted for the state 
and costates in this thesis. 

Although the shape function for the control u only defines a constant value 
within the element, values of u at additional points are available. For instance, 
once the nodal values for the states and costates are found, then one may use the 
optimality condition ( dH/dit = 0) to solve for u at a nodal point. In fact, this is 
how one finds a value for Uf to use in the 8t / equation. Also, if the states and 
costates are approximated by some continuous curve fit through the nodal values 
obtained from the solution, then the control could be approximated at any instant 
in time by using the optimality condition. 
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2.3. Example: A Fixed-Final-Time Problem 


As the first optimal control problem, the transfer of a particle to a rectilinear 
path will be examined (see Fig. 2.2). This is an example taken from article 2.4 of [1]. 
Let and x ( 2 ) denote the position of the particle at a given time and X( 3 ) and X( 4 ) 
denote the particle’s velocity at a given time. (A subscripted number in parentheses 
refers to the state index to avoid confusion with the element index.) The thrust 
angle u is the control and the particle has mass m and a constant acceleration a. 

The state equations are defined as 



rO 

0 

1 

01 


0 1 


0 

0 

0 

1 

X + < 

0 

X — 

0 

0 

0 

0 

a cos u 


.0 

0 

0 

0 . 


< asinu , 


(2.3-1) 


The final time T is fixed and the problem is to maximize the final horizontal 
component of velocity. Thus, 


L = 0 

(l>=[ 0 0 1 OJ Xf 


(2.3-2) 


The optimality condition dH/du = 0 yields an expression for the control of tan u = 
A( 4 )/A( 3 ). There are also two terminal constraints on the states. These are that the 
particle arrive with a fixed final height ( h ) and that the final vertical component 
of velocity be zero. The final horizontal component of position is free. These 
constraints can be stated analytically as 


= 


[0 

i 

0 

o' 

- f h \ 

o 

0 

0 

1 



(2.3-3) 


The initial conditions are x(0) = xo = [0 0 0 0J T . Finally, the unknown 

A/’s are eliminated by writing it in terms of other unknowns. In accordance with 
Eq. (2.1-7) 
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0 


A/= | ^ | (2.3-4) 

The /, Z-, <f > , and boundary conditions to substitute directly into Eq. (2.2-6) 
have now been defined. 

These equations are solved by choosing At; = At = tf/N for all i, expressing 
the Jacobian explicitly and using a Newton-Raphson algorithm. For N = 2, suitable 
initial guesses for the nonlinear iterative procedure can be found by simply choosing 
element values that are not too different from the boundary conditions. The results 
from solving the N — 2 equations are then used to obtain the initial guesses for 
arbitrary N by simple interpolation. In all results obtained to date for this problem, 
no additional steps are necessary to obtain results as accurate as desired. 

Representative numerical results for all four states versus dimensionless time 
t/T are presented in Figs. 2.3 - 2.6. For this example, h = 100, T = 20, and 
4 h/aT 2 = 0.8897018. (This last number is chosen to yield a value of 75° for the 
initial control angle of the exact solution available in [1]) The results for 2, 4, and 
8 elements are plotted against the exact solution. It can easily be seen that N = 8 
gives acceptable results for all the states. Amazingly enough, even the very crude 
2 element mesh yields a decent approximation to the answer. 

In Fig. 2.7, the control angle u versus dimensionless time t/T is presented. Once 
again, the results are seen to be excellent for N = 8. Note the extra data points 
which are available for the control when we make use of the optimality condition 
at the nodal and midpoint values of the elements. This is of great value since it is 
the control variable which is of the most interest. 

Three of the four costates are constants for all time and this method yields two 
of these exactly. The third costate is very close to the exact answer. The fourth 
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costate corresponding to the vertical component of velocity A( 4 ) is shown in Fig. 2.8. 
The results compare nicely with the exact results. 

A plot of the relative error of the performance index J = ^ /,(3) the endpoint 
multiplier v\ versus the number of elements is shown in Fig. 2.9. It is seen to be 
nearly a straight line on a log-log scale. The slope of the line is about —2 which 
indicates that the error varies inversely with the square of iV, similar to a-posieriori 
error bounds as formulated in usual finite element applications [53]. Notice in 
Fig. 2.9 that there is a bend in the endpoint multiplier curve. It is not unusual 
for mixed formulations to have an error curve that is not monotonically decreasing. 
It should be noted that developments of mathematical proofs for convergence and 
expressions for error bounds are not state-of-the-art for mixed methods. However, 
some initial error estimate studies for the weak principle are given in Chapter 9. 
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Horizontal Position/h 


1.5 



Dimensionless Time 


Fig. 2.3: Dimensionless horizontal position x^/h vs. tJT 
Note that the final horizontal component of position is not specified 


27 


Vertical Position/h 



Dimensionless Time 


Fig. 2.4: Dimensionless vertical position x^)/h vs. t/T 
The final height is constrained to be h at the final time 


Horizontal Velocity /(h/T) 



Fig. 2.5: Dimensionless horizontal velocity x^T/h vs. t(T 
Note that the performance index J = £( 3 ) CO 
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Dimensionless Time 


Fig. 2.6: Dimensionless vertical velocity x^T/h vs. t/T 
The final vertical component of velocity is constrained to be zero at the final time 
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Dimensionless Time 


Fig. 2.7: Control angle u vs. t/T 
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Dimensionless Time 


Fig. 2.8: Vertical velocity costate A( 4 ) vs. t/T 
The results for this costate are the least accurate of all the costates 
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Fig. 2.9: Relative error of the performance index 
and the endpoint multiplier vs. N 
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2.4. Example: A Free-Final-Time Problem 


The second optimal control problem is similar to the previous example except 
that now the final time is free and we would like to obtain a given horizontal 
component of velocity ( U ) in the minimum time (see [1], problem 9, article 2.7). 
The algorithm from the preceding example is readily modified to fit this problem 
by noting the following changes. The performance index is now the final time T\ 
so 4> = 0 and L — 1. Also, there is an additional endpoint constraint on the states; 
namely that £(3) = U. With these changes 


and 





(2.4-1) 



'0 10 0 ' 


f M 

^ = 

0 0 10 
0 0 0 1 

Xf - < 

r 
l 0 J 


(2.4-2) 


Along with these changes to the equations, one additional equation is added 
from the coefficient of dtj. The new system of equations is solved in the same 
manner described previously. Again, initial guesses not too far from the boundary 
conditions are satisfactory for N = 2, and these answers are used to obtain initial 
guesses for arbitrary N. 


Representative numerical results for all four states versus dimensionless time 
t/T are presented in Figs. 2.10 - 2.13 for a case with ah/U 2 = 0.75. The results for 
2, 4, and 8 elements are plotted against the exact solution available in [1]. It can 
easily be seen that N = 8 gives acceptable results for all the states. 

The control angle u versus dimensionless time t/T is presented in Fig. 2.14. 
Once again, the results are seen to be excellent for N = 8. In Table 2.1, the initial 
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control u(to) (which can be easily shown to be related to Vi, 1*2, and t/3) and the 
normalized final time ^ are shown to converge quite rapidly as AT is increased. 
Note, however, that the N = 2 and N = 4 approximations for u(t 0 ) are neither 
upper nor lower bounds. This is a common characteristic of mixed formulations. 

In Fig. 2.15, the vertical velocity costate is shown. The agreement of the finite 
element solution and the exact solution is excellent, even for 2 elements. The plot 
of the relative error of the performance index T versus the number of elements N 
is shown in Fig. 2.16. Again, the slope of the line is about —2 indicating that the 
error varies inversely with the square of N. 

Table 2.1: Convergence of u(to) and versus N 


N 

u(fo) (degrees) 

aT 

u 

2 

72.586 

1.8819 

4 

74.736 

1.8531 

8 

75.027 

1.8413 

16 

74.969 

1.8380 

32 

74.950 

1.8372 

exact 

74.944 

1.8369 
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Vertical Position/h 



Fig. 2.11: Dimensionless vertical position x^/h vs. t/T 
The final height is constrained to be h at the final time 
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I. 



Dimensionless Time 


Fig. 2.12: Dimensionless horizontal velocity )/U vs. t/T 
The final horizontal component of velocity is constrained to be U 




Dimensionless Time 


Fig. 2.13: Dimensionless vertical velocity X( 4 )/U vs. t/T 




Thrust 



0.0 0.2 0.4 0.6 


Dimensionless Time 


Fig. 2.14: Control angle u vs. t/T 
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Vertical Velocity Costate 



Dimensionless Time 


Fig. 2.15: Vertical velocity costate A ( 4 ) vs. t/T 
The results for this costate are the least accurate of all the costates 
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CHAPTER 3 


SATURN ONE-STAGE ROCKET MODEL 


In this chapter, a model is presented which is suitable for evaluating the poten- 
tial usefulness of the weak principle for practical problems in optimal control [50]. 
The weak principle will be applied to a one-stage model of the Saturn IB rocket. 
In this case an analytical solution is not available, and the accuracy of the solu- 
tion will be compared with a solution obtained using a multiple shooting method. 
While there would normally exist several inequality constraints in this problem, the 
constraints are not included in this application. However, Chapter 5 contains a 
two-stage Saturn model involving discontinuities in the system equations and in the 
states, and Chapter 8 has a model of an advanced launch vehicle which includes 
staging and inequality constraints. 

3.1. A Four-State Model 

Consider a vehicle confined to vertical plane dynamics and flying over a spher- 
ical, non- rotating earth as depicted in Fig. 3.1. This results in the following state 
model for the states m (mass), h (height), E (energy per unit mass), and 7 (flight- 
path angle): 
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(3.1-1) 


771 

h 

E 

7 


9.81J.p 
V sin 7 
T-D 

m 

T + qSCLc 

mV 


Mt 



cos 7 


where T is the thrust, T vac is the thrust in a vacuum (a constant), D is the drag, 
and V is the velocity. Here a, the angle of attack, has been adopted as a control 
variable. The atmospheric model is given by the following equations: 


p = po ( 1 - 0.00002255/i) 5 ' 256 for h < 11000m 

( h- 11000 \ „ , 

P =Pn exp f — J for h > 11000m 

p ~po exp (— 

a = a 0 Vl - 0.00002255/i for h < 11000m 

a = 295.03ms -2 for h > 11000m 


(3.1-2) 


The aerodynamic and propulsion models are given by the following equations: 


T = T vac - A e p 
r - = R e - \- h 


pV 2 

q ~ 2 

D —qS [C do (M) + a 2 C NQ (M)] 
C La (M)=C Na (M)-C D0 (M) 

M = — 

a 


(3.1-3) 
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The vehicle parameters chosen for this model are based on a Saturn IB launch 
vehicle SA-217 [54] and are 


5= 33.468 m 2 ; = 263.4 s 

(3.1-4) 

T vac = 8155800 N; A e = 8.47 m 2 

The aerodynamic coefficient data Cjv a and Cdq are presented as functions of the 
Mach number M in Tables 3.1 and 3.2 and shown graphically in Figs. 3.2 and 
3.3. The physical constants used in the above model are the earth’s gravitational 
constant p = 3.9906 x 10 14 m 3 s -2 , the earth’s mean radius R e = 6.378 x 10 6 m, 
the sea- level atmospheric pressure po = 101320 Nm -2 , the atmospheric pressure at 
11 km pn = 22637 Nm -2 , the sea-level density of air po = 1.225 kg m~ 3 , and the 
sea-level speed of sound in air do = 340.3 ms -1 . 

The initial conditions specified are m(0) = 5.2 x 10 5 kg, fo(0) = 1800 m, 
£(0) = -6.25 x 10 7 m 2 s -2 , and 7(0) = 75°. The final energy is specified as 
E(t f ) = -4.25 x 10 7 m 2 s -2 . 

The performance index is 


J = $ \ h = m\ tf 


(3.1-5) 


and the final time tf is open. Note that L = 0. 

From Eqs. (2.1-7) and (3.1-5) it is seen that A/ is given by [1 0 u 0J T . 

The quantities necessary for direct substitution into the algebraic equations of 
the weak principle Eq. (2.2-6) have now been defined. 
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Table 3.1: Aerodynamic coefficient C^ a versus Mach number 


M 

C Na 

0.00 

6.20 

0.50 

6.35 

0.98* 

7.70 

1.00 

7.70 

1.02* 

7.70 

2.50 

5.20 

4.40* 

4.70 

5.00 

5.50 

6.00 

6.00 


(* denotes a common end point of two quadratic polynomial curves) 


Table 3.2: Aerodynamic coefficient Coo versus Mach number 


M 

Cd o 

0.20 

1.00 

0.75 

0.45 

0.98* 

0.80 

1.00 

0.80 

1.02* 

0.80 

3.50 

0.20 

6.00 

0.02 


(* denotes a common end point of two quadratic polynomial curves) 
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Mach Number 


Fig. 3.2: Csa versus Mach number 
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Cdo 



Mach Number 


Fig. 3.3: Cdo versus Mach number 
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3.2. Results 


Two codes have been developed that solve the resulting finite element equa- 
tions. One uses the Newton-Raphson method, and the other uses the method of 
Levenberg-Marquardt as coded in the IMSL subroutine ZXSSQ [55]. With the 
former method, the initial conditions need to be reasonably accurate. However, 
running a case for a few elements generates a good approximation for larger num- 
bers of elements as above. Also, with this method one can easily exploit sparsity; the 
computational savings of this will be investigated in later chapters. With the latter 
method, the initial guesses do not need to be very accurate, but the method is not 
nearly as computationally efficient as the former since it generates an approximate 
Jacobian from central differencing. 

In Figs. 3.4 - 3.12, numerical results for the Saturn one-stage model are given. 
In these figures, the finite element results are shown as discrete symbols while the 
solid lines show results obtained from a multiple shooting code, an essentially exact 
solution [56]. In Figs. 3.4 - 3.7, the mass, altitude, specific energy, and flight- 
path angle profiles are shown. For the mass and specific energy, even N = 4 gives 
excellent results. For the altitude and flight-path angle, N = 8 gives good agreement 
with the multiple-shooting results. When higher final energies were tried, the rate 
of change of the energy became very steep, and the finite element results became 
more difficult to obtain. One must remember, however, that these results are not 
realistic because of the absence of state constraints and, particularly for the steep 
energy gradients, the use of single-stage modeling. 

Figs. 3.8 - 3.11 show the costate profiles. The results for the mass costate and 
specific-energy costate axe not quite as accurate as were the results for the states. 
Adding more elements would, however, increase the accuracy of the solution. 

In Fig. 3.12 the angle of attack profile is shown. Here the N = 4 solution agrees 
well with the exact solution except near the beginning of the trajectory. Note the 
smooth and rapid convergence at t = 0 as N increases. 
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Fig. 3.4: Mass versus time 
“exact” and finite element results 


51 




specific energy (m''Z/sec''A> 



Time (sec) 


Fig. 3.6: Specific energy versus time 
“exact” and finite element results 



Fig. 3.7: Flight-path angle versus time 
“exact” and finite element results 
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Mass costate 



Fig. 3.8: Mass costate versus time 
“exact” and finite element results 
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Fig. 3.10: Specific energy costate versus time 
“exact” and finite element results 
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Fig. 3.11: Flight-path angle costate versus time 
“exact” and finite element results 
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CHAPTER 4 


DISCONTINUITIES IN THE SYSTEM EQUATIONS 


Although many practical problems in optimal control theory can be solved by 
using the weak formulation derived in Chapter 2, there are problems that require 
more generality. Specifically, the weak formulation must be extended to allow for 
discontinuities in the states and/or discontinuities in the system equations. Such dis- 
continuities might arise when finding the optimal trajectory of a multistage rocket. 
At the time the first stage is dropped, the mass state would suffer a discontinuity. 
Furthermore, the thrust of the rocket (one of the system parameters) would also 
change at this staging time, thereby creating a discontinuity in the system equa- 
tions. These discontinuities produce jumps in the states, costates, and possibly the 
control variables. 

The derivation of the weak formulation to include state discontinuities (or 
jumps) and discontinuities in the system equations is similar to the derivation in 
Chapter 2, (see also [57]); however, special care must be taken because of the un- 
known staging time. Therefore, the details of the derivation are presented below 
for a problem with one discontinuity. Of course, the extension of the formulation 
to problems with multiple discontinuities is possible and should not cause any con- 
fusion. 
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4.1. Theory for State Discontinuities 


Consider a problem with one discontinuity where the time of discontinuity will 
be called the staging time and denoted by t a . The formulation must be modified to 
accommodate the unknown staging time t a , the constraint on the states at t 3 (as 
opposed to a constraint on the states at the final time), the jump in the states at 
t 3 , the jump in the costates at t s , the change in state equations at t a (due to the 
change in system parameters), and finally the transversality condition to find t 3 . 
Furthermore, the control u may be discontinuous at the staging time. 

Now, let a system be defined by a set of n states x and a set of m controls 
u. Furthermore, let the system be governed by a set of state equations of the form 
x = /j(x,u,£) prior to t 3 and x — /jj(x,u,t) after t 3 . Elements of the performance 
index, J, may be denoted with integrands Lj(a;,u,<) prior to t 3 , Ljj(x,u,t) after 
t s , discrete functions of the states and time <j)[x(t),t] defined at the initial and final 
times to and tf, and discrete functions of the states and time <j> 3 [x(t),t] defined just 
before and after the staging time signified by tj and <+ respectively. In addition, 
any constraints imposed on the states and time at to and t j may be placed in sets 
of functions ^>[x(f),f], whereas constraints imposed at t a and tj" may be placed 
in tl> 3 [x(t),t}. These constraints may be adjoined to the performance index by dis- 
crete Lagrange multipliers v and u 3 respectively. Finally, the state equations may 
be adjoined to the performance index with a set of Lagrange multiplier functions 
A(t) which will be referred to as costates. For notational convenience, we define 
$ = <j>{x(t),t] + v T rp[x(t),t] and = <f> a [x{t),t] + vjip 3 [x(t),t]. Also, a constraint 
equation is adjoined to the performance index (as was done in Chapter 2) to trans- 
form the strong boundary conditions to weak ones. For variable t 3 and tf , this 
yields a performance index of the form 
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(4.1-1) 


J = 



[Lj + A t (/j - x)] dt + 



+ ^ T (/lI _ *)] dt 


+ $s 



* T/ A, \ 

+ a [x — x) 

*0 


tf 

to 


Now, denote with S() the variation of () when holding time fixed, and let d() be 
the variation of () when time is allowed to vary. The fixed and free-time variations 
at t = to and t — tf are related by (see [1] or [2]) 


6x(to) = dx(to) and Sx(tf) — dx(tf) — x\ t} dtf (4.1-2) 


and similarly for A. (Note that to is considered to be a fixed time so that dto = 0.) 
The free and fixed-time variations at t = tj and t — tf are related by 


6x(t s ) = dx{t s ) — x\ t - dt 3 and 6x(t+) = dx(t+) — x\ t +dt 3 (4.1-3) 


and similarly for A. It is noted that if a particular state (or costate) does not have 
a discontinuity at t — t 3 , then the corresponding free-time variation is continuous 
at t 3 , i.e., dx(tj) = dx(tf). Now, proceeding with the development of the weak 
form, the first variation of J is taken and the Sx term appearing in both integrands 
is integrated by parts. The resulting equation is 
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dJ 


= f ’ \sL\ + Sf [ A + S\ T (fi - x) + Sx T a] dt - &r r A| 

J tn 


' <0 
rtf 


T\\t. 

to 


+ f ^6Lji + Sf^X + 6\ T (fii~ x) + Sx T \^ dt-Sx T A||+ 

J tt 


+ dx T 


d$s 

dx 


+ 6ujip s _ + Sv T il> 
t* 

+ Sa T (x — x) * + a T (dx — dx) +dtf 

to to 


+ dx 1 


1 7 



+ 

9 

t f 

to 


(4.1-4) 




L n + A (hi - *) + & 
+ dt 3 1 [Lj + \ T (fi - i)] t _ - [L u + A T (/ n - i)] t + + ^} 


A necessary condition for an extremal of J is that the first variation be zero. 
Also, the admissible variations of the states must be continuous at the initial and 
final times and therefore (dx - dx) \\[ = 0. For notational convenience, define 
H i = L^ + A T /j, ffjj = Xjj + A r /jj, and 


* 1 ,.- 


dx 


and A | 


<0 


dx 


(4.1-5) 


Finally, by using Eqs. (4.1-2) and (4.1-3) to eliminate 6x at t 0 , t 3 , f+, and tf , the 
above equation can be simplified to 
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i: 

/; i 


S\ T (fi — x) + 8x J 


+ UA J (/ n -x) + fo J 


+ ^ujxp s 
+ dx T 


+ bv 1 \j) 
I*/ 


"I - dx 


T 


m 
m 

m T » 


A + 

A + 


+ Su T 

+ Su T 


m 

m 


dt 


dt 


t. 


(^A — Xj * + 6a T (x — x) 




(4.1-6) 


It is easily verified that the necessary conditions for an extremal of J , as defined 
in [1], are contained in the coefficients of the virtual quantities [*.e., the 6() and 
d() terms] of Eq. (4.1-6). This is described in Chapter 2 and thus only the new 
boundary conditions which appear as the result of the staging are discussed here. 
Specifically, the requirement for the coefficient of 8v s to vanish yields Eq. (3.7.3) 
of [1]. The requirement for the coefficient of dx{t ~ ) to vanish yields Eq. (3.7.11), 
whereas the requirement for the coefficient of dx(tf) to vanish yields Eq. (3.7.12). 
Finally the coefficient of dt s is the transversality condition in Eq. (3.7.13). 

Having satisfied the requirement that none of the fundamental equations are 
altered, the weak formulation may now be derived. First, we choose 8a = d\. Next, 
the x and A terms in Eq. (4.1-6) are integrated by parts. Also, Eq. (4.1-2) is used 
to eliminate dx and d\ at to and t j . The resulting equation is 
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C[ 


8\ t /j + 8\ T x — 8 x t A + 6x J 


m 




dt 


+ /!' [« t /h + (^) T + ^ (^r) 


dt 


+ 8uJ %l> 9 


dx 


1*7 


+ 8u t xI>\ 


*/ 

«o 


+ £ar A 


x 


*/ 

to 


+ 6\ T x 


‘ -8x t A 


(£)* 

**'(""**) **-(»i|, ; '»n|,,+ S )- 0 


(4.1-7) 


Note that a (A - A)x term and an (x - x)A term should appear in the dtf equation 
but these terms are each zero in accordance with the natural boundary conditions 
(see the dx and 8a coefficients in Eq. 4.1-6). Eq. (4.1-3) is now used to eliminate 
dx at t~ and <+. Again, the natural boundary conditions in Eq. (4.1-6) must be 
used to avoid changing the dt s coefficient. The weak principle is now given as 


+ 


£ [« r /i + - * r * + (l 1 )" + © 

[’ [« r /ll + ^ + fa* (©" + 6u T ©) 


dt 


dt 


(4.1-8) 


•cl */ 


+ 6x J A -8 A J x 
*0 


<0 


+ 8x 


d$ s 


+ 8u;xp s 


1 7 


+ dt s I H i 


1 7 


dx 

-Hu 


T t+ _ ,*+ 


+ 8 X 1 x 


1 7 


1 7 


+ 8u T rp 


d$ 3 \ , f TT d$\ 

,t + m;) +dtf { Hn + -m) 


= 0 


This is the governing equation for the weak Hamiltonian method for optimal control 
problems of the form specified. It will serve as the basis for the finite element 
discretization described below for constructing of candidate solutions (i.e., solutions 
which satisfy all the necessary conditions). 
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4.2. Finite Element Discretization 


Due to the staging, the finite element discretization must be handled somewhat 
differently than was done in Chapter 2. Therefore, for clarity, full details of the 
discretization will be given below. Let the time interval from to to t~ be broken 
into N\ elements and the time interval from tf to tj be broken into N? elements. 
For notational convenience, define N = N\ + iV 2 . The nodal values of time for these 
elements are t{ for i = 1,2, . . . ,N + 1 where to = ti, is = ijVj+i, and t/ = i- A 
nondimensional elemental time r is defined as 


t -tj _ t-tj 
ti+1 tj Atj 


(4.2-1) 


Since one derivative of Sx and S A appears in Eq. (4.1-8), linear shape functions 
for Sx and 6 A may be chosen. Since no derivatives of x or A appear, then piecewise 
constant shape functions for x and A are chosen. These shape functions are taken 
to be 


and 


Sx = Sxf(l 

~T) + Sx i+l T 

(4.2-2) 


if t = 0; 


X = < Xj 

if 0 < r < 1; 

(4.2-3) 

l *i+l 

if r = 1 



and similarly for 6X and A. The superscripted “+” and signs signify values 
just before and after the subscripted nodal value. For all nodes except the Ni + 1 
node which corresponds to t s , the values for x and A, as well as for 6x and 6A, are 
equal on either side of the node. Furthermore, it is important to understand that 
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£+ = x 0 = x(t 0 ), \~i = A 0 = A (to), x N +i = Xf = x(tf), and X N+1 — A/ — A (</). 
We again choose u — U{ and Su = Sxii as shape functions for the control and its 
variation. 

Plugging in the shape functions described for x, A, and u y substituting t = 
+ rAi;, and carrying out the integration over r from 0 to 1, we obtain a general 
algebraic form of our Hamiltonian weak form for optimal control problems of the 
form specified. The algebraic equations are 


r 

£K 


1=1 


A, + 


AU / dHi 

2 \ dx J i 


-sxf 


— Sx] 


+1 


_ AU (dHi\ 

' 2 (&J, 


T1 


+ sx; 


*+i 


+ 6*T 


- X (hh 

Zi + ^ (/l)i 

H+_ M + 


Atl {id ) . } + (^ +1 " i ^*+ 1 + at.* ) 

-f T 

-8xf\+ + 8\+ T x+ + 8x+ T i+1 {^) + 6X N T ,+i i N T l +i 




~ Sx N l + 1 ( ~^T ) ^A/i+l Z N,+I 


N 


+ £ 

1=^+1 1 

— Sx~~ T 

0X i+l 

+ 8uJ 


A i + 


A,- 


At, f dHji \ T 
2 V dx )i 

au /dH u y 




2 V / f . 

T 


+ 6X 


«+i 


x <-^T (Ai)< 




W*- + £) 


— ^A N+1 X Ar+1 + t ^*l<7 "t ^ V* — 0 


( 4 . 2 - 4 ) 


where H = H(x,u,t) and = H(x,u,t). Note that there are not only boundary 
conditions at to and t/, but also additional boundary conditions at the staging time 
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t s , namely the coefficients of 8 : +1 , <5xjy i+1 , ^jvj+i’ an< ^ ^ Ni+i • This is where 
the jumps in the states and costates are allowed to occur. Thus, rather than the 
nodal values of the shape functions canceling one another at the staging time (as 
they do at all other internal nodes), the nodal values at tj and will be distinct. 

Eq. (4.2-4) is a system of nonlinear algebraic equations. The coefficient of each 
arbitrary virtual quantity ( 8x , 8 A, 8u , dt 3 , dt f , 8 z/j, and 8v ) must be set equal to zero 
in order to satisfy Eq. (4.2-4). However, not all of the virtual quantities above are 
independent. As stated earlier, 8xf — 6x~ for all i except z = Ni + 1, which is the 
node number corresponding to the staging time. At this node, it was observed from 
the calculus of variations that the virtual quantities suffer a discontinuity. Now, the 
coefficients of 6x^ i+1 and <5xjy i+1 may be treated as separate equations; however, 
a different option has been chosen in an attempt to simplify the equations to be 
programmed. Define 


and 


8x N i+1 = ^xjVj+i + 8rj\ 


(4.2-5) 


^iVi+i = 6xn 1+ i - 8r]x (4.2-6) 

When these values are substituted into Eq. (4.2-4), then two new arbitrary 
virtual quantities appear, namely and 8r]\. Fig. 4.1 helps clarify the as- 

sembly process of the virtual states. The figure shows three straight lines depicting 
linear shape functions over three elements. For the non-jump node Ni, the virtu- 
ally quantities are equal at the node and replaced with At the jump node, 

however, <$xjy i+1 and £x^ i+1 are replaced with an average value 8xsi+ 1 - Another 
virtual quantity, 8r)\, also appears but is not shown in Fig. 4.1. The coefficient of 
6x7 ^+ i is now of the same form as all the other 8x terms and the 8tj\ coefficient 
contains an equation to extract the needed nodal value of A at t 9 , namely A^q+i- 
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To simplify matters further though, the coefficient of 8r}\ is replaced with a still 
simpler expression to extract the needed nodal value. This expression comes from 
the following recursive equation derived in Chapter 2. 


Zi+i = 2 Zi - Z{ 


(4.2-7) 


where z represents either the state x or the costate A. The first nodal value z\ is 
equal to the initial values of the states and costates which are represented by xq 
and A 0 in Eq. (4.2-4). The same process is done with the <5A]0 i+1 and <5A)^ i+1 terms 
so that a ^A^-i-i and 8r] x axe introduced. As a final step, all superscripted “+” and 
signs can now be dropped (except on xyvj+i and A^q-fi, because these values 
are still distinct) since they are equal at all nodes but the staging time node which 
has now been handled. The algebraic equations now take the form 
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(4.2-8) 


Eq. (4.2-8) may be used to solve optimal control problems of the form specified. 
Note that for this three-point boundary- value problem the elements must be assem- 
bled over the entire time interval. Only the nodal values (the hatted quantities) of 
the states and costates at to, t 3 , an d t/ appear in the algebraic equations. 

Often in practical applications, these algebraic equations may be simplified 
slightly in order to minimize the number of unknown variables. Consider the case 
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(as will be seen in Chapter 5) where ip 3 is a scalar function at t s and , i,e., there 
is a known jump in one of the states. 

More specifically, since the jump in the states expressed as i^ l+1 — ^Nx+i ls 
known, then the coefficient of + i will be replaced with an actual numerical 
value, such as the drop mass of a booster rocket stage. Also, would be retained 
to define the unknown staging time, but $+ would be set to zero since the jump has 
already been solved for. Then, the jump in the costates (from the natural boundary 
conditions) will be defined as 


A- - A+ - 
a Ni+1 a Ni+1 - Qx 


(4.2-9) 


To use Eq. (4.2-8), one would let x^ i+1 and A^ +1 be the unknown variables 
and then replace +1 (from the physical jump condition) and Ajy i+1 (from Eq. 4.2- 
9) in terms of other unknowns. Note that if there is no jump condition on a 
particular state or costate, then the nodal values just before and after staging are 
equal. Fig. 4.2 helps clarify this process. In Fig. 4.2 are shown three piecewise 
constant shape functions spanning three elements. The solid black circles in the top 
half of the figure represent the hatted or nodal values of the states at the beginning 
and end of the element. Note that only the jump node is labelled. After assembly, 
two conditions can occur. If there is no jump, then the nodal values from the 
left and right sides are equal and do not appear in the algebraic equations. These 
“disappearing” nodal values are represented by the hollow circles in the lower half 
of the figure. (There values are recoverable though as described earlier.) If there is 
a jump, then the nodal values x)y i+1 and are no ^ e 9 ua l> but only £jVi+i * s 

treated as an unknown (represented by the solid black dot) and is eliminated 

from the equations in terms of i)v l +i - Thus it is now one of the “disappearing’ 
nodal values also and depicted by a hollow circle. 
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In addition, one must introduce the optimality condition ( dHjdu = 0) at t~ 
and tj in order to solve for the values of the control u just before and after the 
staging time. Finally, the coefficient of the dt s equation (t.e., the continuity of the 
Hamiltonian) gives the extra equation to solve for the unknown staging time t a . 

Also, as was noted in Chapter 2, for most problems the initial conditions are 
given for all n states and thus, in accordance with Eq. (4.1-5), all the initial costates 
are unknown. Therefore, instead of treating elements of v at t = to as unknowns 
and replacing A| (o with these unknowns, we will instead treat A |* 0 as unknowns and 
eliminate the <5i/| (o equations from the weak principle. 

It may be instructive to count the number of equations and unknowns for a 
given problem. Consider a problem with n states, m controls, qi constraints on the 
states at t a and <72 constraints on the states at tf. There will be JVi elements in the 
first stage and N 2 elements in the second stage. The given boundary conditions will 
be for xo and A /. The number of unknowns is 2n (for Ao and Xf) + 2n(N\ + N 2 ) (for 
X{ and A i in the first and second stages) + 2n (for x^ i+1 and A^ i+1 ) + m(iV 1 + jV 2 ) 
(for iii in the first and second stages) +3m (for u at and if) +q\ (for i / s ) 

+ 92 (for v) + 2 (for t a and t /). The number of equations is 2n(N\ + iV 2 + 1) (for Sxi 
and £Ai) + 2n (for Sr] x and Sr}\) +m(Ni + N 2 ) (for Stii) +3m (for dHjdu = 0 at 
and tf) +91 (for Su s ) +<72 (for & v ) +2 (for dt a and dtf). Thus, the number 
of equations and unknowns is the same. 

Eq. (4.2-8) is actually not much more complicated to program than are the 
equations presented in Chapter 2 [see Eq. (2.2-6)]. In fact, the one-stage rocket 
model was modified for the two-stage rocket model presented in the next chapter 
without any serious complications. One of the most tedious and time-consuming 
tasks was changing the program to account for the linearization of the new and 
unknown staging time, as this variable appears in about half of the equations. 
However, the general code described in Chapter 11 eliminates all the programming 
difficulties. 
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Fig. 4.2: Assembly of states 
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CHAPTER 5 


SATURN TWO-STAGE ROCKET MODEL 


In this chapter, a more complicated and realistic model than that of Chapter 3 
is presented which is suitable for evaluating the potential usefulness of the weak 
Hamiltonian finite element approach for real-time guidance of a launch vehicle. A 
two-stage, four-state vehicle is considered that is simplified by not allowing for any 
inequality constraints. This model allows us to incorporate the theory developed in 
Chapters 2 and 4. 


5.1. The Model 

The same vehicle model from Fig. 3.1 is used. The dynamical equations are 


m = — 


9.81 1 


sp 


h = V sin 7 

e=( t -^ 


7 = 


V 

m J 

T + qSC La 

mV 


(5.1-1) 


“ + (7 - £ v ) COS7 


where T is the thrust, T vac is the thrust in a vacuum, D is the drag, and V is the 
velocity. Here o, the angle of attack, has been adopted as a control variable. 


The atmospheric, aerodynamic, and propulsion models are taken to be the 
same as in Eqs. (3.1-2) and (3.1-3). The vehicle parameters chosen for this model 
are based on a Saturn IB launch vehicle SA-217 [54] and are 


75 


I ap j = 263.4 s; I apn = 430.4 s 

T vaC j = 8155800 N; T vacn = 1186200 N (5.1_2) 

^ej = 8.47 m 2 ; ^ = 5.29 m 2 ; £= 33.468 m 2 

where subscripts “I” and “II” refer to the first and second stages respectively. 

The aerodynamic coefficient data C'y Q and Cdo <ire also taken to be the same as 
the one-stage model. They are presented as functions of the Mach number M in Ta- 
bles 3.1 and 3.2. The graphs are shown in Figs. 3.2 and 3.3. The physical constants 
used in the above model are the earth’s gravitational constant p = 3.9906 x 10 14 
m 3 s 2 , the earth’s mean radius R e = 6.378 x 10 6 m, the sea- level atmospheric pres- 
sure po = 101320 Nm -2 , the atmospheric pressure at 11 km pu = 22637 Nm -2 , 
the sea- level density of air p 0 = 1.225 kg m -3 , and the sea- level speed of sound in 
air ao = 340.3 ms -1 . 

The performance index is 


J = = m| t/ 


(5.1-3) 


and the final time tf is open. The initial conditions specified are m(0) = 5.2 x 10 5 
kg, h{ 0) = 1800 m, E{ 0) = -6.25 x 10 7 m 2 s" 2 , and 7 (0) = 75°. The burnout mass 
of the first stage is 192000 kg ( ift~ = — 192000) and the drop-mass of the 

booster is 51000 kg. The final energy is specified as E(tf ) = —4.25 x 10 7 m 2 s -2 
(V> = £(*/) + 4.25 x 10 7 ). 

From Eqs. (4.1-5) and (5.1-3) it is seen that Ay is given by [1 0 iq 0J T . 
Note that the only jumps are in the mass state and the mass costate, and these 
jumps are 


m(<7) -m(i+) = 51000 

rn{t s ) — ^m(t^) = V s 


(5.1-4) 
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5.2. Results 


The finite element equations are solved using the method of Levenberg- 
Marquardt as coded in the IMSL subroutine ZXSSQ [55]. Running a case for a few 
elements generates a good approximation for larger numbers of elements. Initial 
guesses do not need to be very accurate, but the method is not nearly as compu- 
tationally efficient as a Newton- Raphson procedure where sparsity in the Jacobian 
can be exploited. 

Numerical results for the Saturn two-stage model are given in Figs. 5.1 - 5.10. 
Discrete points are given for 2, 4, and 8 elements in each time interval (denoted by 
N = Ny : N 2 on the graphs). The converged results corresponding to N\ = 8 and 
N 2 = 16 are shown as solid lines. Note that the number of elements in each interval 
is completely arbitrary. 

Figs. 5.1 - 5.4 show the four states. In Fig. 5.1, notice how nicely the jump in 
the mass is allowed for by the discretization. Also, we point out that the awkward 
altitude profile of Fig. 5.2 (*.e., the strange drop at the end of the trajectory) is 
a result of an unrealistic model. The model is unrealistic due to the absence of 
inequality constraints, and due to the large angles of attack (more than 30° at some 
points) even though small angles were assumed in the state equations. However, 
the model does suffice to illustrate the power of the method. 

The four costates are shown in Figs. 5.5 - 5.8. Again the jump is allowed 
for very accurately by the discretization. Also, note that for the N = 2:2 case, 
the jump is actually in the wrong direction. Even though this is a very inaccurate 
result for the mass costate, the N = 2 : 2 case is still close enough to the real 
answer to allow us to interpolate and run higher numbers of elements. This gives 
some indication of the robustness of the method. 

The control (o) is shown in Fig. 5.9. There is a jump in the control at the 
staging time due to the change in the thrust vector magnitude. The jump is solved 
for in the program by enforcing the optimality condition at t~ and tf . Remember 


77 


that the optimality condition is used at all nodes to get the control, and this leads 
to the extra data points on the control profiles. 

As an indication of the accuracy of the method in a global sense, the Hamilto- 
nian was observed to converge to zero (the exact answer) all along the trajectory as 
is seen in Fig. 5.10. The finite element results are converging to the exact solution 
as N increases. 
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Fig. 5.1: Mass versus time 
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Fig. 5.3: Specific energy 
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Fig. 5.5: Mass 
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Fig. 5.7: Specific-energy costate versus time 
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Fig. 5.10: Hamiltonian versus time 
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CHAPTER 6 


CONTROL INEQUALITY CONSTRAINTS 

Many practical optimal control problems have certain constraints imposed on 
the magnitude of the control variables. This is done for a variety of different reasons. 
For instance, if the control is to be produced by a power supply, then there may 
be constraints on the controls so that the power supply does not become saturated. 
Another reason for a control constraint would arise when studying flight vehicles 
where the structural integrity of the vehicle might be jeopardized by too large a 
control variable. 

The weak principle will now be derived to include problems with control in- 
equality constraints. After the derivation is given, a simple example problem is 
presented. The numerical results are compared to the exact solution. Of particu- 
lar interest is the performance in terms of execution time and accuracy versus the 
number of elements used to represent the time span of the problem. 

6.1. General Development 

The same problem statement as was given in Chapter 2 is used for this chapter. 
Namely, consider a system where the states are continuous. Now, suppose that g is 
a p x 1 column matrix of constraints on the controls of the form 

g(x,u,t) < 0 (6-1-1) 
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One way of handling inequality constraints is to use a “slack” variable [58]. 
The idea is that if g < 0 then g plus some positive number (t. e., the slack variable) 
is equal to zero. Thus denoting the slack variable by k 2 , then the following p x 1 
column matrices for K and SK, the variation of K, may be defined. 


K=[k\ k\ ... kl\ T 
SK = \2k\Ski 2&2<5&2 ... 2k p 8kp\ T 


( 6 . 1 - 2 ) 


Now, from Eq. (6.1-1) 


g(x,u,t) + K = 0 


(6.1-3) 


Eq. (6.1-3) will also be adjoined to the performance index J by using p La- 
grange multiplier functions 


<*(0 = Lm /*2 ... lipi T 


(6.1-4) 


The performance index now takes the form: 

f t} 

J = [L(x,u,t) + \ T (f - i) + p T {g + I<)] dt + ^lh +a T (x - x)f/ o (6.1-5) 

■J to 

To derive our weak principle, it is necessary to take the first variation of J and 
set it equal to zero. For simplicity, the derivation below is for the case of fixed final- 
time. The case for free final-time is discussed after this derivation. For notational 
convenience, the following variables are introduced. 
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( 6 . 1 - 6 ) 


. 9 $ 

Ao — 

ox 


and *, = £ 

3 OX 


Also, as is shown in Section 2.1, the Lagrange multiplier a can be chosen so that 
Sa = SX. 

The first variation of J is 



-f 8v T xl>\\[ + ^x t A|^ -f <5A T (x — x)|*' + a T (<5x — Sx) |j' = 0 


(6.1-7) 


The admissible variations to the states must be continuous at the initial and fi- 
nal times and therefore ( Sx — = 0. Furthermore, it is noted that for most 

problems, the initial conditions are given for all n states and thus, in accordance 
with Eq. (6.1-6), all the initial costates are unknown. Therefore, instead of treating 
elements of u at t = to as unknowns and replacing A| to with these unknowns, A |( 0 
will be treated as unknowns and the 5^| to equations will be eliminated from the 
weak principle. Finally, the weak principle is obtained by integrating the i: term 
in Eq. (6.1-7) by parts. Denoting the variations of the variables at the initial and 
final times with subscripts 0 and /, then the resulting equation is 
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( 6 . 1 - 8 ) 



This is the governing equation for the weak Hamiltonian method for fixed-time 
problems with control inequality constraints. It is easily shown by integrating the 
8x and <5 A terms by parts in Eq. (6.1-8) that all the Euler-Lagrange equations are 
the same as in [1] and that all boundary conditions are now of the natural type. 
When the final time is allowed to vary, Eq. (6.1-8) remains unchanged except that 
one term is added, given by 


8t i 




u dt ) 


(6.1-9) 


Note that in this term, values for u are required at t f. To obtain it, one must 
also find the values of K and // at t /. These unknowns are found by setting the 
coefficients of 8uf, 8Kj, and SfiJ equal to zero in the following 



( 6 . 1 - 10 ) 


Thus, the formulation has now been developed to handle fixed and free-time 
problems with inequality constraints on the control. 
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6.2. Finite Element Discretization 


The same finite element discretization described in Chapter 2 is used for this 
problem. The only added step is to define shape functions for A and fi. Since the 
time derivatives of K , // , <5 1 \ , and 8fi do not appear in the formulation, then 


K = I\, n = m 
6K = 8Ki 8n = 8fii 


( 6 . 2 - 1 ) 


By substituting Eq. (4.2-1) and the shape functions described above into 
Eq. (6.1-8), and carrying out the element quadrature over r from 0 to 1, a gen- 
eral algebraic form of the Hamiltonian weak principle is obtained. Again, if the 
time t does not appear explicitly in the problem formulation then all integration 
is exact and can be done by inspection. If t does appear explicitly, then t may be 
approximated by a constant value over each element and the integration may still 
be done by inspection. The latter case occurs in the example problem presented 
shortly. For N elements, there are 2 n(N + 1) + mN + q + 2 Np equations and 
2 n(iV + 2) + mN + q + 2 Np unknowns. Therefore, 2n of the 4n endpoint values for 
the states and costates (x 0 ,A 0 ,Xf, and A/) must be specified. In general, x 0 (the 
initial conditions) is known in accordance with physical constraints. Also, A / can 
be specified in terms of other unknowns with the use of Eq. (6.1-6). Now there are 
the same number of equations as unknowns. These equations may be used for any 
optimal control problem of the form specified. 

Although the shape functions for u, K , and ft only define a constant value 
within the element, values for these variables at additional points are also available. 
For instance, once the nodal values for the states and costates are found, then one 
may use the optimality condition ( dHjdu = 0), the constraint equations (g + A = 
0), and the condition that either k or n be zero (k/x = 0) to solve for u, I<, and g, 
at a nodal point. This procedure is used in the following example problem. 
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6.3. Example 


This example is taken from section 3.8 of [1], The problem is to minimize 


J = 



(6.3-1) 


where T = 10, x and u are scalars, and the initial condition is x 0 = —19.945596. 
The state equation is 


x = h(t)u with h(t ) = 1 + t — — t 2 


(6.3-2) 


The following two control inequality constraints are imposed. 


gi = u — 1 < 0 

92 = ~(u + 1 ) < 0 

The exact solution is found to be x(T) = -17/39 and for the control 


(6.3-3) 


u(t) = 


-x(T)h(t) for 0 < t < 2 
1 for 2 <t < 11/3 

— x(T)h(t ) for 11/3 < t < 8 
-1 for 8 < t < 10 


(6.3-4) 


The algebraic equations which come from the weak principle can be verified to 
be 


A 0 - A (1) = 0 

A (i ) - A (<+1) = 0 for i = 1, 2, . . . , JV — 1 (6.3-5) 

A^> -x f = 0 


for the Sx coefficients, 
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for * = 1,2, ...,1V (6.3-6) 


S«> + A (»h po] + $ _ $ = 0 ' 

fciV, 0 = 0 

k^fi^ = 0 ► 
« (i > - 1 + k[ i)7 = 0 
-«(*■> - 1 + 4 <)S = 0 , 


for the 6u, SK , and coefficients, and 


x< i+ D - - — 

2 


{a [*«] 

-x< N) 


u (i) + h [^ i+1) ] u ( ‘ +1) } = 0 

"T*^] ^ N) + i f = ° 


for 


i — 1,2,.. 


N- 1 


(6.3-7) 


for the SX coefficients. Note that is an average time value for the i th element 
and (if A ti = At = tf/N for all i ) can be expressed as 


fi» = — — ^At for * = 1,2, ... ,1V 


(6.3-8) 


Recall from [1] that one of the additional necessary conditions for problems 
with control constraints is that the multipliers be greater than or equal to zero for 
a minimizing problem. Therefore, in practice, the multipliers // appearing in the 
first of Eqs. (6.3-6) are squared to ensure their positivity. Further recall that if the 
constraint is not violated, then fi = 0. This condition is satisfied by the second 
and third equation in Eq. (6.3-6) which implies that either k or // is zero for each 
element. 
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It is readily apparent from the equations shown in Eq. (6.3-5) that all the 
costate variables are equal to xj. Therefore, these equations were eliminated and 
all costates that occurred in the remaining equations were replaced with Xf. The 
remaining 6N + 1 algebraic equations were solved using a Newton-Raphson method 
and a FORTRAN code written on a SUN 3/260. The sparse, linearized equations 
are solved using subroutine MA28 from the Harwell subroutine library [59]. This 
subroutine takes advantage of sparsity which leads to great computational savings. 

Table 6.1 shows the convergence rate of xj — x(T), the elapsed computer 
time for the first 5 iterations, and the percentage of zeroes in the Jacobian (be., 
the sparsity) versus the number of elements. The x(T) column shows that the 32 
element case has almost converged on the exact solution. Note further that the 
approximate ,r(T) is not an upper bound of the exact value, which is common in 
mixed formulations. The third column of Table 6.1 gives the elapsed computer time 
for five iterations. It is easily seen that there is a modest increase in computer time 
with an increase in the number of elements. Note that in some cases a converged 
answer is found in five or fewer iterations. This is because the answers obtained 
from a small number of elements (say 2 or 4) may be interpolated to generate initial 
guesses for a higher number of elements. Thus, it is possible to solve a 16 or 32 
element case in about 1.5 seconds. Finally, the extremely sparse structure of the 
Jacobian is demonstrated in the last column. This strongly encourages the use of 
a “smart” sparse matrix solver such as MA28. This subroutine leads to quicker 
solutions and tremendous savings in memory allocation since only the nonzeroes of 
the Jacobian need be stored. 
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Table 6.1: x(T), elapsed computer time and percent sparsity of 
Jacobian versus the number of elements N 


N 

x(T) 

Time (sec) 

Sparsity (%) 

1 

-4.0632 

0.42 

65.3 

2 

-.82795 

0.44 

80.5 

4 

-.44065 

0.66 

89.6 

8 

-.43360 

0.76 

94.6 

16 

-.43928 

1.03 

97.3 

32 

-.43588 

1.52 

98.6 

Exact 

-.43590 

— 

— 


Results for the control u are shown in Fig. 6.1 for 2, 4, and 8 elements and 
the exact solution. Note that although the 2 element case does not define the 
constraint boundaries very accurately, it is accurate enough to generate guesses for 
the 4 element case. Thus, in a problem with many constrained and unconstrained 
arcs, a small number of elements could still be used to generate guesses for a higher 
number of elements. Also, it is interesting to note that as few as 4 elements have 
essentially converged on the exact solution. 
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CHAPTER 7 


STATE INEQUALITY CONSTRAINTS 

Optimal control problems with state inequality constraints are in general very 
difficult to solve. Although the problem does not seem more difficult conceptually 
than problems with control inequality constraints, there are additional necessary 
conditions to reckon with. These additional necessary conditions are a result of what 
are sometimes referred to as “tangency” conditions [1]. These tangency conditions 
arise from the following physical considerations. 

If the constraints are of the form S(x,t ) < 0, then successive total time deriva- 
tives of S are taken and f(x,u,t ) is substituted for x until an expression explicitly 
dependent on the control u is obtained. If p total time derivatives are required, 
then S is called a pth-order state variable inequality constraint. Now, since S(x,t) 
can be controlled only by changing its pth time derivative, no finite control will 
keep the system on the constraint boundary if the path entering onto the constraint 
boundary does not meet the “tangency” conditions. These conditions are that S 
and all the time derivatives of S up to p — 1 are zero. These conditions also apply 
to the path leaving the constraint boundary. 

There have been several papers over the past 30 years presenting necessary 
conditions of optimality for optimal control problems with state-variable inequality 
constraints. In [60], the authors adjoined the pth derivative of the constraint to the 
performance index and allowed the tangency conditions stated above to form a set 
of interior boundary conditions. These boundary conditions require discontinuities 
in the costates at the junction points between constrained and unconstrained arcs. 
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However, one may arbitrarily pick the entry point as the place to satisfy these 
boundary conditions, and therefore the costates and Hamiltonian are discontinuous 
at the entry point and continuous at the exit point. 

A new set of necessary conditions was given in [61]. Therein, the constraint 
was adjoined directly to the performance index. It was shown that unconstrained 
arcs had to satisfy tangency constraints at both ends of a constrained arc. 

Solution of problems with state constraints may be handled in a variety of 
ways. For example, a penalty function approach was presented in [62], A Valentine 
transformation technique (as was used in Chapter 6) was presented in [63]. This idea 
transformed a constrained problem into an equivalent unconstrained problem by 
adding additional state equations. Also, the control variable is transformed and the 
new control appears linearly in the formulation. The weak principle cannot handle 
linearly appearing controls, so this idea was not used. Finally, [64] demonstrates 
how a state constraint present in the full-order problem may be transformed to a 
control constraint in the reduced-order problem. 

New necessary conditions for optimal control problems with state inequality 
constraints were developed in [65]. Therein, the authors tactfully say: “We do not 
imply that the necessary conditions obtained by previous workers are incorrect, but 
rather, that, inasmuch as they underspecify the conditions at the junction, there 
exists the possibility of non-stationary solutions satisfying these conditions . . ..” The 
authors generated an admittedly contrived example where the necessary conditions 
of [l], which are identical to those in [60], were satisfied by a non-extremal solution. 
These new necessary conditions are summarized below. 

The new necessary conditions derived in [65] may be summarized as follows. 
Rather than take successive time derivatives of the constraint S(x,t ) until the con- 
trol appears explicitly, the authors adjoin S directly to the Hamiltonian, as was 
done with control inequality constraints in Chapter 6. The Hamiltonian now takes 
the form 
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H =L + \ T f + r)S(x,t) 


(7-1) 


where 77 , S, and u are scalars. 

As with control inequality constraints, the multiplier rj is positive if S' = 0 and 
zero if S < 0. At junction points t{ of boundary and interior arcs, the costates may 
be discontinuous. The boundary conditions are 

= A((-) - (If ) t ; K<i)>0 (7-2) 

and, in addition, 

H(tf) = H(t-)-v si - (7-3) 

An extremely interesting consequence of the new conditions is pointed out in 
the paper for problems which possess a Hamiltonian which is said to be regular. 
The Hamiltonian H is regular if along a given trajectory, H has a unique minimum. 
In this case, it was shown that the control u and its ( p — 2) time derivatives are 
all continuous. Now, the interesting consequence of the new conditions is that for 
an odd-order constraint greater than two, the trajectory will, at most, only touch 
the boundary if the ( p — l)th derivative of u is discontinuous at the junction point. 
Note that for p = 1 , the control is continuous, so that boundary arcs are permitted 
for the first-order case. 
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7.1. General Development 


Consider once again a system as defined in Chapter 2. Now, suppose that there 
is a pth order scalar constraint on the states and time defined by S(x,t ) < 0. The 
first attempt to apply the present methodology to problems with state inequality 
constraints made use of the necessary conditions presented in [65]. These necessary 
conditions lead to successful and accurate solution strategies for states that only 
touch ( i.e ., do not ride) the constraint boundary. As is derived in [65], for constraints 
of odd order greater than one, the solution can at most only touch the constraint 
boundary if the Hamiltonian is regular. However, for cases where the states ride 
the constraint boundaries for a nonzero length of time, the algebraic equations 
developed by the weak form are singular. Private discussions with Jason Speyer and 
Dan Moerder indicate that the cause is related to a reduced-dimensional manifold; 
however, we have not been able to develop a nonsingular weak form as of now. 

Below are presented two very similar weak formulations using the necessary 
conditions of [65] for touch-point cases and [l] for ride cases. Fortunately, the nec- 
essary conditions presented in [l] are accurate for first and second order constraints 
where the solution often rides the constraint boundary and the conditions in [65] 
are accurate for the touch-point cases. It is noted that most practical applications 
will be third-order or less. 

7.1.1. Touch-Point Cases 

The weak formulation is now derived for touch-point cases. Assume that there 
is only one touch-point over the time interval of interest whose time will be denoted 
by t ip . In this case, the state constraint is nothing more than an interior boundary 
point which creates a jump in the costate. 

The performance index J now takes the form: 
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(7. 1.1-1) 


J — / [L(x,u,t) + X T (f — i)] dt + I [L(x,u,t) + A T (/ — i)] dt 

J <0 

+ + $l!' + oc T (x - *)|Jj 

To derive the weak principle, it is necessary to take the first variation of J and 
set it equal to zero. This variation, and the entire development of the weak principle 
is almost identical to the derivation given in Chapter 4. The only difference is that 
the state equations are the same on either side of the constraint. As usual, we 
introduce 


Ao 



and 



(7.1. 1-2) 


Also, as is shown in Chapter 4, the Lagrange multiplier a can be chosen so that 
8 a = dX. The final form of the weak principle is obtained after integrating by 
parts so that no derivatives of the states or costates appear. After defining the 
Hamiltonian H = L + A T f and denoting the variations of the variables at the 
initial, touch-point, and final times with subscripts 0, 1, and / respectively, then 
the resulting equation is 


+ 


-Sx T X + SX T f + 6X T x + 6x 
-Sx T X + SX T f + SX T x + 8x 




rttj 

Jt 0 

[— Sx 1 X + 8X 1 f + 8X 1 x + Sx 1 + 8u 

+ Sv T s \u r +6uT $\t, +6x T v \ 


Ts , c\Tr , c \ TV , fSHY , t .. r 


( 

m 


dt 


dt 


+ 8x T fXj — <5xjfA 0 — SXjxf + 8 XqX 0 

dS 


-f- dt 


tp 


H(tT p ) - H(t+) + 


+ dt f 




di 

a7J 


(7.1. 1-3) 
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This is the governing equation for the weak Hamiltonian method for problems with 
touch-point state inequality constraints. It is easily shown by integrating the 8x 
and <5 A terms by parts in Eq. (7. 1.1-3) that all the Euler-Lagrange equations are 
the same as in [65] and that all boundary conditions are now of the natural type. 

One simplification may be made to Eq. (7. 1.1-3). If the control is continuous 
across f< p (as is guaranteed if the Hamiltonian is regular), then it is possible to 
simplify the dt tp equation since then f{tj p ) = f(tf p ) = f(tt P ) and L{t^ p ) = L(tf p ) = 
L(t tp ). From the necessary conditions that are found in [65] or from the ones that 
could be found from Eq. (7. 1.1-3), it is seen that 


A T (f tp ) - A T (tt p ) = 


Tu+ 


dS 

dx 


(7.1. 1-4) 


Now, rewriting the coefficient of dtt p as 


H(t7 p ) - H(tf p ) + V 1 


dS 

dt 


[A r (*< p ) - A r (tJ,)] /(< tp ) + vi 


dS 

dt 


dS 


dS 


dS 


(7. 1.1-5) 




we see that the condition for continuity of the Hamiltonian reduces to the condition 
that the first total time derivative of the constraint be zero at t tp if the control is 
continuous. 

An example of a touch-point case is given in the next section. 
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7.1.2. Boundary Arc Case 


For cases where there is a boundary arc (i.e., the solution rides the constraint 
boundary for a nonzero length of time), then the weak formulation must be modi- 
fied. For simplicity, consider the case where the solution has an unconstrained arc 
between to and < en , followed by a constrained arc between t en and t ex , and then 
another unconstrained arc between t ex and tf. Introducing a new Lagrange multi- 
plier function 77 to adjoin the pth derivative of the constraint S to the performance 
index, then J becomes 


J= / [L(x,u,t) + \ T (f-i)] 

Jto 

+ It 


dt 


+ A T (/ -X) + T} 


dPS~ 

Up 


dt 


(7. 1.2-1) 


+ J ’ [L(x, u, t ) + A T (/ - i)] dt + UiN\t en + + ot T {x - x)| 


*0 


where N is a column matrix defined as 


f - £3] (71-2-2) 

Analogous steps to those described in Chapter 4 lead to a weak formulation for 
state constraint problems which ride the constraint boundary. There are only two 
minor differences. One is that the time line must be discretized between t 0 and t en , 
from t en to tex, and from t ex to tf. Also, there will be equations corresponding 
to the 6r) coefficient over the interval from t en to t ex , just as there were with the 
control constraint of Chapter 6. These equations are that the pth derivative of S 
be zero between t en and t ex . 

Also note that the dt en and dt ex equation will reduce to the condition that 
the pth total time derivative of 5 be zero. The proof is identical to that shown in 
Eq. (7.1. 1-5) after noting that 
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(7. 1.2-3) 


A T (C») - ^(4) = V, 


dN 

dx 


Two examples of problems with a boundary axe are given in the next section. 


7.2. Example: A First-Order Problem 

Consider the classical brachistochrone problem in Section 3.11 of [1]. Let x 
and y define the horizontal and vertical (positive downward) position of the particle 
respectively. The governing state equations are 


i = (2 9V)' n COS u 
y = (2jy) 1/2 sinu 


(7.2-1) 


where g is the acceleration due to gravity (a constant) and the control u is the 
angle that the tangent to the path makes with the horizontal. The problem is to 
minimize the time it takes for the particle to move from the origin to any point on 
the line x = L. The optimal path to the unconstrained problem was found over 200 
years ago by several mathematicians to be a cycloid. The problem takes a new twist 
though when a state inequality constraint is added. Let S = y — x tan 9 — h < 0 
where 8 and h are constants. The first total time derivative of S yields 


S = (2 gy) 1 ! 2 sin(u — 9)/ cos 9 = 0 (7.2-2) 

or u = 9 along the constraint boundary. Thus, this is a first order constraint. 

By adjoining S to the performance index, the solution was readily found, al- 
though the programming was a little tricky because of the three unknown times to 
reckon with. Fig. 7.1 shows the trajectory of the particle and Fig. 7.2 shows the 
control history. In the figures, (2:2:2) designates that 2 elements were used on each 
of the unconstrained arcs and the constrained arc in between. The finite element 
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solution for 2 and 4 elements in each phase are graphed versus the exact solution 
found in [1]. The answers are excellent for the states, but are a little off on the 
controls. When more elements are run, the answer is seen to converge on the exact 
answer. Fig. 7.3 shows the log of the error in the entry, exit and final times versus 
the number of elements. The lines are all approximately straight with a slope of 
about —1. This indicates that the times have an error proportional to about 1/N. 
This is uncharacteristically inaccurate as compared to other results presented thus 
fax. 
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Fig. 7.2: it versus time 
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Fig. 7.3: Log of the error versus the number of elements 


110 


7.3. Example: A Second-Order Problem 


This example is taken from section 3.11 of [1]. The problem is to minimize 


J =\fr' dt 

The state equations are 

x i = u 

±2 = Xi 


(7.3-1) 


(7.3-2) 


The state inequality constraint S(x,t) = x 2 — t < 0 is to be imposed. For certain 
values of £, the solution only touches the boundary, whereas for other values of t 
the solution rides the boundary. 

The algebraic equations for both examples were solved using a Newton-Raphson 
method and a FORTRAN code written on a SUN 3/260. The sparse, linearized 
equations are solved using subroutine MA28 from the Harwell subroutine library 
[62]. 

The state x 2 is shown in Fig. 7.4 for the single touch-point case. Results for 
2, 4, and 8 elements on either side of the touch-point (denoted by 2:2, etc.) are 
compared to the exact solution. Note that even the 2:2 element case lies essentially 
on the exact solution. In Fig. 7.5, the state x 2 is shown for an example case where 
the state rides the boundary. Here, there are three time intervals and the number 
of elements in each interval is denoted by 2:2:2 etc. Again we see that the 2.2.2 case 
has essentially converged on the exact solution. 

One drawback of the weak formulation is that two separate codes had to be 
written to solve this problem. Also, one must determine in advance if the solution 
will ride or just touch the constraint. However, with the general code described in 
Chapter 11, this is a simple and quick thing to do. 
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Fig. 7.4: x 2 versus time for a touch-point case of £ = 0.2 
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Time (secs) 


Fig. 7.5: X 2 versus time for a boundary arc case of £ — 0.1 
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CHAPTER 8 


AN ADVANCED LAUNCH VEHICLE 


As stated in Chapter 1, future space transportation and deployment needs are 
critically dependent on the development of reliable and economical launch vehicles 
that will provide flexible, routine access to orbit. The objective of the Advanced 
Launch System program is to place large payloads - 100,000 to 150,000 pounds - 
into low Earth orbit at an order of magnitude lower cost per pound. The program 
also seeks to make the entire scope of space launch operations significantly more 
routine compared to present methods and procedures that are highly dependent on 
ground operations. The goals of the ALS program can only be met by development 
of reliable and efficient on-board algorithms that are capable of calculating real-time 
optimal trajectories. 

In this chapter, a model of an advanced launch vehicle is presented [66]. This 
is a two-stage, four-state vehicle with control and state inequality constraints to 
be imposed. The results from the finite element algorithm are seen to compare 
favorably with multiple shooting results. Of great interest is the fact that the 
answers are obtained in a very short time interval as compared to the time span of 
the entire trajectory. 
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8.1. A Model for an Advanced Launch Vehicle 


A two-stage, four-state vehicle is considered that has two control inequality 
constraints and one state constraint. (Only one of these constraints is violated and 
is therefore the only one included in the results.) 

We confine our attention to vertical plane dynamics of a vehicle flying over a 
spherical, non-rotating earth (see Fig. 3.1). This results in the following model for 
the states m (mass), h (height), V (velocity), and 7 (flight-path angle): 


m = — 

Qisp 

h = V sin 7 
^ T cos a — D 

7 = 


H sin 7 


( 8 . 1 - 1 ) 


m 

' T sin a + L 
mV 


M’-*) 


cos 7 


where T is the thrust, T vac is the thrust in a vacuum, D is the drag, and L is the 
lift. Here a, the angle of attack, has been adopted as a control variable. 


Note that now V is being used as a state instead of E. There are two reasons 
for this. First, using V as a state simplified the algebraic equations. Secondly, since 
E weis two or three orders of magnitude larger than the other states in the Saturn 
models, convergence was easier to obtain with V as a state. 

The aerodynamic and propulsion models are given by the following equations: 


T = T vac - A e p(h ); M = 


V 

a(h) 


r =R e + h] q = 


pV 2 


D = qSC d(M, a); L = qSC L (M , a) 


( 8 . 1 - 2 ) 


The atmospheric data for density, pressure, and speed of sound are obtained 
from the 1975 standard atmospheric data [67]. 
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The vehicle parameters chosen for this model are 


Si = 131.34m 2 ; % = 65.67 m 2 

T vac = 25813400 N; T vac = 7744020 N 

1 11 (8.1-3) 

■A e j = 37.515 m 2 ; ^ e n = H-254m 2 

I 3p j = 430.0 s 

where subscripts “I” and “II” refer to the first and second stages respectively. 

The aerodynamic coefficient data C d and Cl are functions of the Mach number 
M and angle of attack a. The physical constants used in the above model are the 
earth’s gravitational constant fi = 3.9906 x 10 14 m 3 s -2 , the earth’s mean radius 
R e = 6.378 X io 6 m and the acceleration due to gravity g = 9.81 ms 2 . 

The three constraints are 

gi(x,u,t ) — aq — 2925 rad-Pa < 0 

g 2 (x,u,t)= — (aq + 2925) rad-Pa < 0 (8.1-4) 

q(h, V) — 40698.2 Pa < 0 

where only the first constraint will be enforced since the other two are not violated. 
The performance index is 


J = $ \ tf = m\ t] (8.1-5) 

and the final time tf is open. The initial conditions specified are m(0) = 1.52345 x 
10 6 kg, h(0) = 400 m, V(0) = 64.48941 m/s, and 7(0) = 89.5°. The final conditions 
are h(t f ) = 148160.0 m, V(t f ) = 7858.1995 m/s, and j(t f ) = 0.0°. The burnout 
mass of the first stage is 645500 kg and the drop-mass of the booster is 98880 kg. 
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8.2. Computational Aspects 


The code for this model was written on a SUN 3/260. The code was written 
to be as efficient as possible so as to get a feel for the actual run-times one might 
see in an on-board computational setting. 

An explicit Jacobian is formed within the code and a Newton- Raphson method 
was employed. The Jacobian was 85% to 95% sparse for the runs made. Great 
computational savings came from taking advantage of the Harwell sparse matrix 
solver MA28AD [59]. The code was run with double precision. 

Initial guesses were, of course, necessary for the Newton-Raphson method. A 
Taylor series approach was taken that generated initial guesses for all variables. 
These guesses, although crude, were good enough to make a converged run with 
boundary conditions that differed somewhat from the specified conditions. We 
were then able to slowly “move” the boundary conditions out to the specified ones. 
An explanation of the Taylor series approach is given in Chapter 10. 

8.3. Results 


In Figs. 8.1 — 8.8, numerical results for the ALV model with no constraints 
enforced are given for 2, 4, and 8 elements per time interval, where the number 
of elements is denoted by ( A'i : jVj ) on the plots. These results are compared 
to a multiple- shooting code as a check on the accuracy of the method and of the 
program. The four states are shown in Figs. 8.1 - 8.4 and the costates are shown 
in Figs. 8.5 - 8.8. For all cases, the (8:8) run lies on the essentially exact curve 
corresponding to the multiple shooting (MS) code. In general, even the (4:4) run 
yields an excellent approximation to the solution. 

The control is shown in Fig. 8.9. Although the (8:8) run is close to the exact 
curve, it has not converged on the answer. Due to the large slopes and sharp peaks 
in the control, the finite element method required 24 elements in the first stage to 
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converge on the solution. However, it is important to note that we were still able 
to run only 8 elements in the second stage. Thus, it may be feasible to cluster the 
elements to refine the solution. 

Fig. 8.10 shows a graph of aq. (The dynamic pressure q is not shown because 
it does not violate the constraint.) It can be seen clearly that only one control 
constraint, < 71 , is violated. This is the only constraint that was added to the program 
to generate the next two graphs. 

The (4:4) results above were obtained in 5.5 CPU seconds on a SUN 3/260, 
and five iterations were required with a Newton-Raphson method. Of course, the 
number of iterations depends on the quality of the initial guesses. In an on-board 
computational setting, the initial guesses should be pretty good since they would 
probably be determined from a previously obtained solution. 

For Figs. 8.11 and 8 . 12 , the control constraint g\ was included into the computer 
model. Since the constraint is just barely violated, there was essentially no change in 
the data when the multiple shooting code was run; therefore, it is not seen in these 
two graphs. However, for illustrative purposes, the unconstrained case, the realistic 
constraint, and two unrealistic constraints are shown for the finite element case. 
Even for the lowest of the constraints, the states, costates, and dynamic pressure 
are virtually unchanged. Also, no significant extra computer time was expended. 

Finally, as a feel for the global convergence of the method, the Hamiltonian of 
the unconstrained system is plotted in Fig. 8.13 for 2, 4, 8 , and 16 elements per 
stage. There is a nice convergence toward the exact answer of zero with an increase 
in the number of elements. 
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Mass (kg) 



Fig. 8.1: Mass versus time 
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Altitude (m) 



Fig. 8.2: Altitude versus time 


120 





Fig. 8.3: Velocity versus time 
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Fig. 8.4: Fligl 
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Mass Costate 



Fig. 8.5: Mass costate versus time 
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Fig. 8.6: Altitude costate versus time 
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Fig. 8.7: Velocity costate versus time 
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Fig. 8.8: Flight-path angle costate versus time 
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CHAPTER 9 


ERROR ESTIMATES 


Numerical examples have shown that the finite element method presented in 
this thesis yields very accurate solutions to initial-value ordinary differential equa- 
tions, dynamics problems, and optimal control problems. It is the intent of this 
chapter to find a relationship between the step size At and the error of the integra- 
tion performed [68] . Error estimates for dynamics have been performed by Hodges 
and Hou [69]. Therein, the authors find that the error for the linear oscillator prob- 
lem to be It is noted that to obtain general error estimates for mixed methods 
is not state of the art. 

The first part of this chapter will be concerned with comparing the local trun- 
cation error of Euler’s method, a second-order Runge-Kutta method, and the finite 
element method for the initial- value ordinary differential equation x = f(x, t). The 
approximation to x(t) will first be given for each of the three methods, and the 
Taylor Series expansion of each approximation will be derived. Then, the local 
truncation error of the three methods will be compared to the exact solution for 
four different /(x,t)’s. The finite element method will be seen to be proportional 
to At 3 for rill examples. 

The second part of the chapter will involve studying the error of the finite 
element method in optimal control problems. A control problem will be solved 
using the weak form developed in Chapter 2. It is possible to develop equations for 
all the unknown variables in terms of At. Consequently, convergence to the exact 
solution is proven and the error in each variable is shown to be proportional to At 2 . 
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9.1. Initial- Value Problems 


Consider the differential equation 


x = f(x,t ) 


(9.1-1) 


where a: is a scalar, t is the time, and x 0 is the given initial condition. It is desired 
to find the local truncation error involved with the integration of this differential 
equation. Several approximation methods are explored. These are Euler’s method 
(also known as a first-order Taylor Series method), a second-order Runge-Kutta 
method, and the finite element method which is developed in Appendix B. Following 
the summary of these methods, four example problems will be examined. 

9.1.1. Methods of Solution 

Let A t be a small time step and let the value of x at t = A t be denoted by x. To 
compare the error of an approximation method with the exact answer, x will be 
expanded in a Taylor Series in At about At = 0. 

The value of x using Euler’s method [70] will be denoted by xem and is defined 
as 


zem = ro + /(fo,0)At 


(9. 1.1-1) 


This method is already in a Taylor Series form where only first-order terms in At 
are present. 

The value of x using a second-order Runge-Kutta method [70] will be denoted 
by xrk • Defining g = f(x 0 , 0), then ^rk is given as 
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(9.1. 1-2) 


ZRK = + At/ 



At 

2 


) 


For clarity, define 


/rk = / 




(9.1. 1-3) 


where it is noted that /rk is always an explicit function of At. The Taylor Series 
expansion of Eq. (9. 1.1-2) is 


£rk « xo + /rk(£o>0)A< + 2 


5/rk(^o,0) At 2 _5 2 /rk(^o,0) At 3 


dAt 


+ 3 


dAt 2 


(9.1. 1-4) 


An expression for the approximation of x using the finite element method, xfe, 
is not quite as easy to produce. This is because the finite element method is an 
implicit method. The equations derived in Appendix B (and Chapter 2 also, but 
for a different reason) are 


x - — f(x,t ) = x 0 
xfe + xo = 2x 

Eq. (9. 1.1-5) may be rewritten as one equation involving only nodal (i.e., hatted) 
values. This equation is 


(9.1. 1-5) 


..sfXFE + io At 
x FE ~ A tf( , — ) = Xo 


(9.1. 1-6) 
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In order to find a Taylor Series approximation for xfe in the above equation, it is 
necessary to rewrite xfe as a polynomial in A t with unknown coefficients. Thus, 
let 


iFE = 5>'+ lA(i (9. 1.1-7) 

«=0 


Eq. (9. 1.1-7) is now substituted into Eq. (9. 1.1-6) resulting in 


G = f>+. Ai' - A tf |^) - x„ = 0 (9. 1.1-8) 

i=0 ' ' 


Again for clarity, define 


/fe = / 


xo + 



(9. 1.1-9) 


where /fe is an explicit function of At. To find the unknown coefficients x \ , x 2 , 
etc., the Taylor Series expansion of Eq. (9. 1.1-8) will be taken about At = 0, and 
the coefficients of each At term will be set equal to zero. The required derivatives 
of G with respect to At, denoted by superscripted numbers, are 
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oo 

G = ^x,+iA t l - Atf F E - xq 

i=0 

OO 

G (1) = Y, ixi+lAV - 1 - /fe - A tf™ 

1=1 

OO 

G (2) = £*(* - 1)^+1 A^" 2 - 2/pe - AtfW 

i=2 


(9.1.1-10) 


G (n) = ^n!x l+1 At*- n - n/pg _1) - A</<£ 

i=n 

Since G' and the derivatives of G when evaluated at At — 0 are the coefficients in 
the Taylor Series approximation of G', then the unknown coefficients for xfe are 
found by setting each equation above equal to zero. This results in 

G L=o = ° = *i -*o 
G ^^Iai=o = 0 = ^2 — /fe|a<=o 

G(2) Ia<=0 = 0 = 2x 3 - 2/pg |a<=o (9.1.1-11) 


& n) \ AM = 0 = n!«„+i - (n - 

Solving Eq. (9.1.1-11) for the unknown coefficients results in 


Xi = Xo 

^2 = /feU(=o 


£ Tl 


r(n-2), 

/ FE I 


(n-2)! JFE 


A (=0 


(9.1.1-12) 
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So, finally, the approximation to x using the finite element method is given by 
Eq. (9. 1.1-7) where the coefficients of the polynomial axe given by Eq. (9.1.1-12). 

The Taylor Series approximation of the Euler method, a second-order Runge- 
Kutta method, and the finite element method have now been derived. These equa- 
tions axe grouped together below for convenience. The coefficients in Eq. (9.1.1-12) 
have been substituted into Eq. (9. 1.1-7). 

xem = xo + f(x o ,0)At 

XRK ~ x 0 + /RK|At=oAt + f^\At=0^i 2 + 2-^RK l A t=oAf 3 (9.1.1-13) 

XFE ~ Xo + /fe|a<= 0^ + /pE |a<= 0 At^ + 2 ^FeI a< = 0 ^ 3 

where 

/RK = /(*» + ^/(io.O),^) (9.1.1-14) 

and 

£q + Xj+iA^_ At 

2 ’ 2 


/fe = / 




(9.1.1-15) 
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9.1.2. Example Problems 


Four example problems with exact solutions will now be examined to compare 
the local truncation error of the three methods described above. It will be shown 
that the finite element method is as good or better than the second-order Runge- 
Kutta method for all the example problems. 

The first example will be a simple linear differential equation in x. Suppose 
= x and xq = 1 so that Eq. (9.1-1) becomes 


x = x 


(9.1. 2-1) 


The exact answer is x(f) — exp(t) so that xex> the exact value of x at t = A t is 
given as 


A/ 3 

xex = exp(Af) « 1 T A< -f 1 


(9. 1.2-2) 


By using Eqs. (9.1.1-14) and (9.1.1-15), f RK and / FE are found to be 


/rk = / 
/fe = / 




2 



1 + iq-n At* 
2 


(9. 1.2-3) 


and the unknown coefficients above are found from Eq. (9.1.1-12). Now, Eq. (9. IT- 
13) may be used to find the Taylor Series representation for each method. These 
expressions are 
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£em = 1 + At 


£rk = 1 + At H — — (9.1. 2-4) 

A t 2 At 3 

ipE ~ 1 + At H — — — I — — 

The local truncation error of the three integration schemes is found by compar- 
ing terms of the Taylor Series in Eq. (9.1.2— 4) with the exact solution in Eq. (9.1.2— 
2). The errors, denoted by e, are found to be 

for Euler’s method 

for a Runge-Kutta method (9. 1.2-5) 

for the finite element method 


e EM 

e RK 

e FE 


At 2 

2 

At 3 

IT 

At 3 

~12 


Thus, it is easily seen that the finite element method is a superior integration scheme 
to the Euler method or second-order Runge-Kutta method for this linear differential 
equation. 

The second example involves a nonlinear differential equation in x only. Let 


x = 1 + cos x 


(9. 1.2-6) 


with xo = 0. The exact solution to this equation is 

x(t) = 2 tan -1 1 


(9. 1.2-7) 
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Following the same procedure as in the first example of finding /rk, /fEi and the 
unknown x’s, then Eq. (9.1.1-13) and the Taylor Series of Eq. (9. 1.2-7) may be 
used to write down the following expressions. 


zem = 2At 

xrk ss 2 At — — At 3 

z 

t-fe & 2 At — —A t 3 

£ 

2 , 

xex ^ 2 At — —At 3 

O 


The errors are found to be 


2At 3 

e EM = — ~ — 

A< 3 

e R K = — 

At 3 

e FE = ~7T“ 


(9. 1.2-8) 


(9. 1.2-9) 


For this case, the finite element method and the second-order Runge-Kutta method 
have the same error, and the error is proportional to At 3 . 

The third example deals with a function of x and t. Consider 


x 


X 2 cos t 


with x o = 1 


(9.1.2-10) 


The exact solution for this problem is x(t) = (1 — sint) 1 . The Taylor Series 
approximations for the three approximation methods and the exact solution axe 
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XEM = 1 + At 


irk ~ 1 + At + A f 2 + 
ife « 1 + At 4- At 2 + 
xex ~ 1 4 A f + At 2 + 


At 3 

8 

9Af 3 

8 

5Af 3 

6 


(9.1.2-11) 


and the errors are found to be 


CEM = At 2 

17 At 3 

e RK = ~~ 24 ~~ (9.1.2-12) 

7At 3 


Once again, the finite element method is more accurate than the second-order 
Runge-Kutta method. The local truncation error is on the order of At 3 . 

For the fourth and final example, consider the special case of Eq. (9.1-1) where 
f is an integrable function of t only. Then, Eq. (9.1—1) becomes 


i = f(t) 


(9.1.2-13) 


Eq. (9. 1.1-6) for the finite element method becomes 

* ,, At 
£fe = £o + A</(— ) 


(9.1.2-14) 


Note that this is now an explicit equation for ife a-nd the Taylor Series expansion 
may be found without specifying a particular /(f). The Taylor Series expansion for 
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Euler’s method, the Runge-Kutta method, and the finite element method are found 
from Eqs. (9.1.1—13) and (9.1.2—14). Note that /rk = /fe = /(Ai/2) so that 

/rk = /fe =-^/ <n, (Y) (9.1.2-15) 

where again the superscripts refer to derivatives with respect to At. 

These expressions, along with the expansion of the exact solution, and the local 
truncation error associated with each method are given in Table 9.1. A summary 
of the errors for each example problem studied in this section is given in Table 9.2. 
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Table 9.1: Taylor Series and error versus exact for x = f(t) 


Method Taylor Series Expansion 


Error 


Euler xo + f(0)At 

Runge-Kutta x 0 + /(0)Ai + A g ' 

finite element x 0 + /(0)At + ^£~ 

Exact x 0 + /(0)A< + ® ^ 


9/(0) At 2 
9At 2 
3 2 /(0) At 3 
9At 2 24 
9 2 /(0) A< 3 
9At 2 24 


Table 9.2: Summary of errors 


/(*><) 

e EM 

e RK 

eFE 

X 

At 2 / 2 

At 3 /6 

At 3 f\2 

1 + COS X 

2At 3 /3 

At 3 /6 

At 3 / 6 

X 2 COS / 

At 2 

17At 3 /24 

7At 3 /24 

/(*) 

/( 1 >(0)A/ 2 /2 

/< 2 )(0)At 3 /24 

/( 2 >(0)At 3 /24 
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The weak Hamiltonian finite element formulation is a mixed method, which is to 
say that there is more than one field variable. Developments of mathematical proofs 
of convergence and expressions for error bounds are not state-of-the-art for mixed 
methods; however, some feel for the accuracy of the weak Hamiltonian formulation 
may be obtained by studying a simple example problem. Afterwards, an error plot 
from the fixed-time problem of Chapter 2 will be examined. 

Consider the following optimal control problem where x and u are scalars. 


J = \ x 0-f + J \ u2 dt 

x — tu with ;r(0) = 4 


(9.2-1) 


The exact solution is readily found to be 

A (<) = 3 

x (t) =4 -t 3 (9.2-2) 

u(t) — —3 1 


It is desired to find an error estimate of each unknown variable using the finite 
element formulation. Since the error estimate will be a function of the element size 
At, then the number of elements N = 1/At will be left as a free parameter. 

The equations to be solved are easily verified to be 
A 0 - A( J) = 0 

A (i) — A Ci+1) = 0 for i = 1,2,...,N -1 (9.2-3) 

\W-x f = 0 
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for the costates, 


u(*> + A(*)^') = 0 for i = l,2,...,N (9.2-4) 


for the controls, and 


*(») - — t* 1 ^ = 4 

2 


^(i+l) _ S (0 _ — ^ 
2 


*««<■'> + ^+I) s (i+1)] = o for i = 1, 2, . . . , N - 1 


At 


.j(»).W fi W + i /S :0 


(9.2-5) 


for the states. Note that is an average time value for the element and can 
be expressed as 


= for i = i >2 ,...,JV (9.2-6) 

2 

There are 3 N + 2 equations with 3 N + 2 unknowns which are and 

«(*) for i = 1,2 ,...,iV, and xj and A 0 . As will be shown, it is possible to find 
an expression for x j solely as a function of At. Then, all the unknowns can be 
expressed as a function of At and the errors of each unknown can be found as 
compared to the exact solution. 

To start with, the N + 1 equations in Eq. (9.2-3) will be solved for the costates 
yielding 


An = A (i) = 


A(*> = x 


Xf 


(9.2-7) 
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Next, the control variables are found from Eq. (9.2-4) to be 


u (0 = = -XftS^ for i = 1,2,. . . ,N (9.2-8) 


After substituting Eq. (9.2-8) for and Eq. (9.2-6) for £**) into Eq. (9.2-5), then 
what remains are N + 1 equations for the N + 1 unknown states. These equations 
axe 


x^ + if f = 4 
/AA 3 

*(H-U _ 5 (0 + i / ) (8i 2 + 2) = 0 for i = 1 , 2 ,..., _ i (9.2-9) 

/AA 3 

-x (iV) +x/f— J (2N - l) 2 +x f =0 


Solving the middle of the above equations for x^ in terms of if results in 
/A t\ 3N ~ x 

* (,) = x w + x f J ^(8; 2 + 2) = 0 for i = 1,2 , N - 1 (9.2-10) 


Using Eq. (9.2-10) to find an expression for it is now possible to solve the first 
and last of Eq. (9.2-9) for x^ and £/. The equation for x/ is 


x f = 


(A!) 3 [ (2 JV - 1) 2 + 1 + EfeT‘(8i 2 + 2)] + 1 


(9.2-11) 
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This messy expression is of more value than it may seem at first. By a series of 
simplifications, it will be possible to prove that if converges to the exact solution 
as the number of elements TV approaches infinity. Furthermore, it will be possible 
to find an estimate of the error for each of the unknowns. 

The simplification will begin with the denominator of Eq. (9.2-11). Noting 
that At = 1/7V, the denominator becomes 


-V 

2 N J 


N - 1 


(27V — l) 2 + 1 + J2(Si 2 +2) 


1=1 


+ 1 


(9.2-12) 


For the case of TV > 1, Eq. (9.2-12) becomes 



N - 1 


47V + 2 + 2(7V - 1) + 8 z’ 2 ) + 1 


i=l 


(9.2-13) 


Making use of the following identity 


N - 1 



(TV - 1)N[2(N - 1) + 1] 

6 


in Eq. (9.2-13), then simplification results in 


1 

8?P 


s -n>- 2 -n) + 1 = 1 


1 

127V 2 


(9.2-14) 


(9.2-15) 


Now, after replacing TV = 1/At in the above equation, then Eq. (9.2-11) becomes 


x f = 


4 

4 

3 12 


(9.2-16) 
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Note that as N — * oo, At — » 0 so Xf = 3 which is the exact answer. Thus, at 
least for this problem, convergence is guaranteed as more and more elements are 
taken. 

To study the error in xj , the Taylor Series expansion of Eq. (9.2-16) in At will 
be taken about At = 0. The resulting equation is 


Xf 


o 3 At 2 
+ 8 ~ 2 ~ 


9 At 3 

16 1T 


(9.2-17) 


The exact value of xj is 3, so the error is At 2 . The final value of the state is the 
least accurate of all the state variables. 

Since all the costate variables are the same and equal to xj (see Eq. 9.2-7), 
then the error in the costates for this problem is equal to 3At 2 /16. Finally, consider 
the error in the control variable at the point t = 0.5. This is a non-moving node 
(for N even) and allows for an error estimate of the control. From Eqs. (9.2-8) and 
(9.2-17), 


u = —Xf/2 


3 A , 9 , 

-1.5- —At 2 + At 3 

32 192 


(9.2-18) 


The exact control is uex = — 1.5 , so the error is 

3A t 2 

IT < 9 - 2 - 19 > 

In summary, convergence of this simple control problem has been proven as the 
number of elements increases. In addition, it has been shown that the error in the 
costates and the final value of the state is proportioned to At 2 . The error for the 
control is also proportional to At 2 . 
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It would be difficult, if not impossible, to use this procedure on a more com- 
plicated problem; however, from observed numerical results, it is believed that the 
finite element method is second-order (At 2 ) accurate for almost all problems. To 
help support this statement, consider once again the fixed-time trajectory problem 
presented in Chapter 2. Fig. 2.8 (repeated in Fig. 9.1) shows a plot of the relative 
error of the performance index J = x and the endpoint multiplier V\ versus the 
number of elements. The relative error of the final control value is also included. 
In Chapter 2, it was noted that the slope of the line is about —2 which indicates 
that the error varies inversely with the square of N , or that the error is propor- 
tional to At 2 . This is similar to a-posteriori error bounds as formulated in usual 
finite-element applications [53]. 
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Number of Elements 


Fig. 9.1: Relative error of the performance index, i/i, and final control vs. N 
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CHAPTER 10 


INITIAL GUESSES 


Almost all of the problems solved to date have used a Newton-Raphson proce- 
dure to solve the nonlinear algebraic equations. This procedure, of course, requires 
initial guesses. And not all initial guesses will converge to the solution. 

Several different ideas have been used throughout the course of this research. 
For example, in the two examples presented in Chapter 2, initial guesses were chosen 
that were not too different from the boundary conditions given in the problem. 
This method worked well for these simple problems and converged solutions were 
obtained without much difficulty. Of course, this idea is nothing more than a trial 
and error method. 

As the problems studied became more difficult, initial guesses became harder to 
obtain. This was due in part to the fact that the costates have little physical meaning 
and hence their range of possible values are generally unknown. For the Saturn 
one-stage model, the solution was obtained by using the method of Levenberg- 
Marquardt [55] which did not require very good guesses. The two-stage solution was 
obtained by “slowly” perturbing the vehicle parameters and boundary conditions 
and resolving the problem at every perturbation. Solutions for the actual specified 
conditions were obtained fairly rapidly using this procedure. 
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10.1. Taylor Series Approach 


With the attempt to automate the program to handle different vehicle pa- 
rameters and boundary conditions more easily, a new initial guess procedure was 
developed. This approach is based on the Taylor series expansion of one of the 
states. The method is outlined below. 

Recall from Eq. (3.1-1) of the Saturn one-stage model the following equations 
for m (mass), h (altitude), E (specific energy), and 7 (flight-path angle): 

m =f\{h) 
h=f 2 (h,E, 7 ) 

( 10 . 1 - 1 ) 

E = f 3 (m,h,E,j,u) 

7 = h, E, 7 , u) 


The first element discretized equations for the states are 


m(1) - =™o 

S (1, -f/1 ,, =A 0 

E^-fff>=Eo 

7 (1) -Y/1‘ ) =7o 


( 10 . 1 - 2 ) 


where aK 1 ) denotes the midpoint value of the state x in the first element, At is an 
arbitrary time step, / (1) is the value of / evaluated at the midpoint of the first 
element, and fo is the given initial condition of the state x. 

Using a Taylor series expansion for h (U results in 


« ho + —f2(ho,Eo,-yo) 


(10.1-3) 
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Now an actual value for h ^ has been obtained. If Eq. (10.1 3) is substituted 
into the second equation of Eq. (10.1-2), then the resulting equation after simplifi- 
cation is 


7 (1) ) = h(ho,Eo,7o) (10.1-4) 

The first, third and fourth equations from Eq. (10.1-2), along with the above 
equation are four equations with the four unknowns (ra^\ E^\ and u^). 

Once the midpoint values are found, the nodal values are extracted and the process 
is repeated. We thus time march until the boundary conditions are approximately 
satisfied. Initial guesses on the costates are then obtained by time marching back- 
wards from the final boundary conditions. The only unknown boundary condition 
for the costates was found by evaluating the Hamiltonian at the final time. 

The above procedure did generate guesses for all the unknowns and the guesses 
did indeed converge on the solution. 

Unfortunately, there are several serious drawbacks to this idea. First, there is 
no way to control the boundary conditions that are to be satisfied. In other words, 
even if the final conditions on the states were changed, the same initial guesses 
would still be generated. Second, the optimality condition is not satisfied. The 
control solved for is not optimal and is, in fact, nothing more than a parameter to 
satisfy the dynamical equations. Third, and most important, there is no guarantee 
whatsoever that the initial guesses will converge. Since a fully automated code is 
desired, a different method to generate initial guesses must be found. 
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10.2. Homotopy and Continuation Methods 

An algorithm for computing Brouwer fixed points wherein the homotopy pro- 
cedure has guaranteed convergence is proposed in [71]. The homotopy method 
was used to generate a numerical algorithm in [72], Therein, the author develops 
“beautifully simple” analytical steps that allow one to find roots of equations. 

The method involves following the zero curve of the homotopy map 
Pa( A, x) = A [x - f(x)] + (1 - A)(x - a) 

starting from (0,a). The zero curve is parameterized by arc length s and it can be 
shown that the zero curve of p a emanating from (0,a) is the solution of an initial 
value problem. When the solution of the initial value problem reaches A = 1, the 
corresponding x is a fixed point of /. Virtually any a will lead to a root of /. The 
art of this procedure lies in following the zero curve accurately, but not too closely 
as this leads to increased computational expense. 

The homotopy method has been applied to nonlinear two-point boundary value 
problems [73,74]. Also, [75] proposes a homotopy algorithm for sparse systems 
of nonlinear equations, which is precisely what the weak principle yields. The 
algorithm leads to substantial computational savings. 

Homotopy methods are rather involved algorithms. When dealing with well- 
behaved equations that have unique solutions, simpler methods may often be used to 
find fixed points of nonlinear algebraic equations. Therefore, a continuation method 
[76] has been adopted to generate initial guesses. The method is now described. 

The most general set of n simultaneous algebraic equations in n unknowns 
Xi , . . . , x n can be written as 


/,-(xx,...,x n ) = 0 for i = l,...,n (10.2-1) 
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Let yi(r), . . . , y„(r) be a set of n functions of a variable r , with 0 < r < 1 and take 


yi (Q) = k, for i = (10.2 2) 


where k { is a constant and selected arbitrarily. Now, require that yi(r), . . . , y n (r) 
satisfy the equations 


£(yi,-..,yn) = /i(*i,...,fcn)(i-' r ) for i = i,. 


. n 


(10.2-3) 


Then, as the right-hand side of Eq. (10.2-3) vanishes at r = 1, the functions 
yi(r),...,y„(r) satisfy, at r = 1, precisely the same equations as do n, . . . ,x„ 
[see Eq. (10.2-1)]. Now, yi(l), . . . , y»(l), and, hence x u ...,x n , may be found 
as follows: Differentiate Eq. (10.2-3) with respect to r, thus obtaining the set of 
first-order differential equations 


df t dy\ dfi dyn 

dyi dr dy n dr 




dfn dyi dfn dyn 

dyi dr dy n dr 


-/ n (fci,...,fc„) 


(10.2-4) 


and perform a numerical integration of these equations using Eq. (10.2 2) as initial 
conditions and terminating the integration at r = 1. 

As was stated previously, ki,...,k n may be assigned any values whatsoever. 
However, it can occur that, for certain choices of h , . . . , k n , some of yi , . . . , y n do not 
possess real values for some values of r in the interval 0 < r < 1, in which event the 
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numerical integration of the differential equations cannot be carried to completion. 
When this happens, one simply changes one or more of . . . , k n . In general, 
results are obtained most expeditiously when k \ , . . . , k n are good approximations to 
Xi, . . . ,x„, respectively. Fortunately, in connection with physical problems, one can 
often make good guesses regarding xi,...,x n , and hence assign suitable values to 
^ 1 , • • • , k n . Finally, it is worth noting that many distinct sets of values of Jfcj, . . . , Jfc n 
can lead to the same values of x% , . . . , x n . 

Let’s look at a simple example of the use of the above procedure. Consider the 
scalar equation 


f(x) = x — cos x = 0 


(10.2-5) 


Now, let y(0) — k — 0 so that Eq. (10.2-3) takes the form 


f(y) = y - cos y = f(k)( 1 - r) = T - 1 (10.2-6) 

The equation corresponding to Eq. (10.2-4) is 

(l+siny)^ = l (10.2-7) 

so that 

^ = (l+sin y)~ l (10.2-8) 


A second-order Runge-Kutta method is used to numerically integrate Eq. (10.2- 
8). A time step of 0.02 was used. When the integration was completed at r = 1, 
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then y(l), and thus x, was found to have a value of 0.7390654, which is correct to 
the fifth decimal place. 

This continuation method is used in the general code described in the next 
chapter. In practice, the integration is performed and then one or two Newton- 
Raphson iterations are required to obtain the final answers. This could be avoided 
by using a more accurate integration scheme, but this tends to lead to increased 
computational effort. 
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CHAPTER 11 


AN ALGORITHM FOR OPTIMAL CONTROL PROBLEMS 


The weak principle for optimal control problems has been fully developed in 
Chapters 2, 4, 6, and 7. The formulation is capable of solving optimal control 
problems that have continuous states, costates, and controls, and problems with 
discontinuities arising from staging (i.e., discontinuities in the system equations), 
control inequality constraints and state inequality constraints. The algebraic equa- 
tions which come from the weak formulation may be derived prior to specifying 
the problem to be solved. It is this feature in particular that allows for a general 
problem-solving environment to be created. 

The main goal of the general code is to reliably solve a large class of optimal 
control problems with a minimum of user interaction. Specifically, it is desired to 
create an environment where the user does not have to write subroutines. To this 
end, a general code has been developed on a SUN 3/260 workstation and requires a 
FORTRAN 77 compiler, MACSYMA [77], and the Harwell subroutine library [59]. 
The general procedure can be broken into three parts that must interface together. 
The first part is the FORTRAN code. This code contains all the subroutines nec- 
essary to solve any of the optimal control problems described above. However, if 
certain problems require table look-up routines (such as aerodynamic data for a 
rocket model), then these subroutines must be given by the user and interfaced to 
the rest of the general code. Thus, there may be a need for some user programming 
for certain problems. The second part of the general procedure is the use of MAC- 
SYMA. The user must supply an input file specifying the problem. This input file 
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is in symbolic form and will be loaded into MACSYMA. MACSYMA will then eval- 
uate all the necessary expressions and automatically generate the FORTRAN code. 
This code is spliced into a template file and becomes one of the subroutines. The 
third and final part of the general procedure will consists of subroutines to generate 
initial guesses that will reliably converge. The continuation method described in 
Chapter 10 is being used. This method converts the algebraic equations to initial- 
value ordinary differential equations. A second-order Runge-Kutta method is used 
to integrate the equations and obtain initial guesses for a Newton-Raphson method. 
This method has worked on all the problems tested to date. 

Every example problem in this thesis has been solved using the general code. 
Although this makes the code useful in and of its self, we wish to emphasize that the 
setup time required to solve these problems is the important factor. Setup time is the 
time required for a user to write all the proper input files and subroutines in order 
to run the program. The setup time to use any of the existing codes (see Chapter 1 
for a discussion of these) can range from several hours to weeks. However, the setup 
time using the general code is only about 10 minutes for any problem (assuming 
that the subroutines for table look-up data are already available). This is because 
only a symbolic input file is required. MACSYMA and the FORTRAN subroutines 
do all the rest. 

As mentioned above, the code can handle a wide range of problems. There 
are some limitations though. One limitation is that the code can currently only 
handle a problem with one or two stages. Also, the code can only handle a scalar 
state constraint. These restrictions are present because of the extra programming 
involved to remove them coupled with the fact that there is no way of testing the 
results analytically. It is felt that there are possibly better ways to handle state 
constraints by using a canonical form as is done in [39]. This is a good problem for 
future research (see Chapter 12). 
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Several example input files are now given to demonstrate the use of the general 
code. Afterwards, a sample output is given for a third-order state constraint prob- 
lem. The chapter concludes with a chart of the times taken to solve the example 
problems. 


11.1. Example Input Files 

The last four pages of this section contain four example input files. These files 
contain four distinct cases to be discussed now. 

Consider the free time trajectory optimization problem presented in Section 2.4. 
The first input file is used to solve this problem. The user is required to supply the 
number of states NS, the number of control constraints NP (zero in this example), 
the number of phases NPH (to be described shortly), the number of controls M, 
and the number of constraints on the states at the final time Q. The next series of 
lines from F[l] to F[4] define the system equations as given in Eqs. (2.3-1). After 
the equations are formed, the user supplies the performance index L and PHI. The 
variable S contains the state constraint, which is zero for this example since there 
are no constraints to be imposed. Then the Q constraints are given in PSI and the 
initial conditions are given in IC. Next the user supplies the final time TF and a 
guess at the value of the final time TFGUES. Since the final time is unknown, TF 
is set to zero and the user gives a guess at the final time. This guess only needs 
to be within an order of magnitude generally. Also, guesses for the states at the 
midpoint of the trajectory and the final point axe given in XGUES. These guesses 
may be very crude and can even be zero for many problems. Since the final value 
of three of the states were known for this problem, crude guesses were easily and 
obviously obtained. Finally, the number of elements to be run is given in NE. 

Regardless of the value of NE, the code automatically starts with the two 
element case and uses the continuation method of [76] and the Newton-Raphson 
method to solve the problem. The code then interpolates the solution to this case 
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and runs a four element case using only the Newton- Raphson method. The code 
continues in this manner until NE is met. If the Newton- Raphson fails to converge 
for the four or higher element case (which is rare) then the program will start that 
case over and try the continuation method to solve the four element case. 

The output was verified to be identical to the solution previously obtained. 
Again, we emphasize that not only is the method accurate as has been observed 
throughout this work, but now the general code can produce these results with a 
simple 30 line input file. This file was created in about five minutes and, as will be 
seen in the last section, answers for 2, 4, and 8 elements were obtained about 6.25 
minutes thereafter. 

The second input file is used to solve the control constraint example of Chap- 
ter 6. Most of the input is the same as in the preceding example. There are, 
however, a couple of important changes to make. One is that now the number of 
control constraints, NP, is 2. Also, these constraints are put in the F array at the 
end of the state equations. Thus, we see that F[2] and F[3] contain the two control 
constraints. Also note that this problem has an explicit dependence on time which 
is handled by the code. Finally, since this problem is a fixed time problem, TF is 
set to the known value and TFGUES is also set to this value. These values must 
be equal for the code to know it is a fixed time problem. 

The third input file introduces another new problem that the code can handle. 
This file will solve the advanced launch vehicle problem presented in Chapter 8. 
For simplicity only in describing this file, the atmospheric effects have been ne- 
glected since there are no analytical expressions for atmospheric and aerodynamic 
conditions. This example, although somewhat unrealistic due to the absence of at- 
mospheric effects, does demonstrate the power and versatility of the general code. 
Several new variables are introduced that need to be explained. First, 1ST appears 
in the sixth line and is set to a value of 1 if there is a second stage involved. The 
next group of lines (7 — 16) are values defined only to simplify the writing of the 
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state equations. F[l] - F [4] define the state equations in the first stage and Fl[l] - 
FI [4] define the state equations in the second stage. The staging time is determined 
by PSITS, which is a constraint on the states at the staging time. Also, the known 
jump in the state is given in JUMP. The index number of JUMP tells the code 
which state is experiencing the discontinuity. Notice again the very crude guesses 
given for the states. It was essential that the code work with a minimum of physical 
insight into the problem. 

As a final example file, a third-order state constraint problem is set up. This 
example is found in [65] and contains an analytical solution. A finite element code 
had not been previously written to solve this problem. Fortunately, the general 
code was able to solve the problem in a matter of minutes. 

A couple of new features are present in this input file. One is that the number 
of phases NPH is something other than 1 for a change. The number of phases is used 
to indicate how often the trajectory enters or leaves a state constraint boundary. 
For this problem, the solution touches the boundary once, so NPH equals 2. Also 
new is an expression for S. S contains the state constraint equation to be enforced 
by the program. It is also found (from results not shown) that a lower (or stricter) 
state constraint causes two touch-points, in which case the user should choose NPH 
= 3. The user is given the task of determining the number of phases in a problem. If 
the incorrect number of phases are chosen that either the constraint will be violated 
or the multipliers which are supposed to be positive will be negative. It is hoped 
that this may be avoided eventually as the general code is further developed. The 
output of this example is given in the next section. 
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NS : 4 ; 

NP : 0 ; 

NPH : 1 ; 

M: 1; 

Q : 3 ; 

F [1] : 1 . 12397 *COS (U (1) ) ; 
F [2] : 1. 12397*SIN(U(1) ) ; 
F [3] : X ( 1 ) ; 

F [4 ] : X (2) ; 

L : 1 . 0 ; 

S : 0 ; 

PHI: 0.0; 

PSI[1] :X(1) -12 . 2129; 

PSI [2] : X (2) ; 

PSI [3] : X (4) -100 . 0 ; 

IC[1] :0.0; 

IC [2] :0 .0; 

IC [3] :0 .0; 

IC [ 4 ] : 0 . 0 ; 

TF : 0 . 0 ; 

TFGUES : 10 . 0 ; 

XGUES [1,1] : 6 . 0 ; 

XGUES [1,2] :12 .2129; 
XGUES [2,1] : 1 . 0 ; 

XGUES [2,2] : 0 . 0 ; 

XGUES [3,1] : 50 . 0; 

XGUES [3,2] : 100.0; 

XGUES [4,1] : 50 . 0 ; 

XGUES [4,2] : 100.0; 

NE : 8 ; 
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NS: 1; 

NP : 2 ; 

NPH : 1 ; 

M: 1; 

Q : 0 ; 

F[l] : (1+T-3*T**2/17) *U(1) ; 
F [2] :0 (1) -1; 

F[3] : -U ( 1 ) - 1 ; 

S : 0 . 0 ; 

L:0.5*U(1)**2; 

PHI : 0 . 5*X (1 ) **2 ; 

IC[1] :-19. 945596; 

TF : 1 0 . 0 ; 

TFGUES : 10 . 0 ; 

XGUES [1,1] : -10 . 0; 

XGUES [1,2] ; -1 . 0 ; 

NE : 8 ; 
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NS : 4 ; 

NP : 0 ; 

M: 1; 

Q:3; 

NPH : 1 ; 

1ST : 1 ; 

TVAC1 : 2 . 58134E07; 

TVAC2 : 7 . 74402E06; 

AE1 : 37 . 515; 

AE2 : 11 . 254; 

SI: 131. 34; 

S2 : 65 . 67; 

ISP ; 430 . 0 ; 

GRAV : 9 . 8 1 ; 

MU : 3 . 9906E14 ; 

RE: 63 78000.0; 

H (X) := RE+X (2 ) ; 

F [1] : -TVAC1/ (GRAV* ISP) ; 

F (2 ] :X (3) *SIN (X (4 ) ) ; 

F [3] :TVAC1*C0S (U(l) ) /X (1) - MU*SIN (X (4 ) ) /H (X) A 2 ; 

F [ 4 ] :TVAC1*SIN (U(l) )/ (X(1)*X(3) ) 

+ (X (3 ) /H (X) -MU/ (X (3 ) *H (X) **2) ) *COS(X(4) ) ; 

Fl [ 1 ] :-TVAC2/ (GRAV*ISP) ; 

FI [2] :X (3) *SIN (X ( 4 ) ) ; 

Fl ( 3 ] : TVAC2 *COS (U ( 1 ) ) /X ( 1 ) - MU*SIN (X (4 ) ) /H (X) A 2; 
Fl [ 4 ] :TVAC2*SIN (U(l) ) / (X(l) *X(3) ) 

+ (X(3) /H (X) -MU/ (X(3) *H (X) **2) ) *COS (X(4) ) ; 

L : 0 . 0 ; 

PHI :X (1) ; 

S : 0 . 0 ; 

PSITS :X(1)- 64 5500.0; 

JUMP [ 1 ] : 9 8 8 8 0 . 0 ; 

PSI [1] :X(2) -148160; 

PSI [2] :X(3) -7858 . 1995; 

PSI [ 3 ) : X ( 4 ) ; 

IC[1] : 1 . 52345E06; 

IC [2] : 4 0 0 . 0 ; 

IC [3] : 160 . 0 ; 

IC [4] : 1 . 4 ; 

TF : 0 . 0 ; 

TFGUES : 300.0; 

XGUES [ 1 / 1 ] : 0 . 7E06; 

XGUES [1,2] : 10000 . 0; 

XGUES [2,1] : 70000. 0; 

XGUES [ 2 , 2 ] : 14 8160 . 0 ; 

XGUES [3,1] : 3000.0 ; 

XGUES [3,2] : 7858 .1995; 

XGUES [4,1] : 0 . 5 ; 

XGUES [ 4 , 2 ] : . 0 ; 

NE : 8 ; 
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NS : 3 ; 

NP : 0 ; 

NPH : 2 ; 

M: 1; 

Q: 3; 

F[l] : X (2) ; 

F ( 2 ] : X (3) ; 

F [3 ] : U ( 1 ) ; 

S :X (1) -0.3; 

L:0. 5*0(1) **2 ; 
PHI: 0.0; 

PSI [1] :X(1) ; 

PSI [2] : X (2) +1; 
PSI [3] :X(3) -2; 

IC [ 1 ] : 0 . 0 ; 

IC [ 2 ] : 1 . 0 ; 

IC [3] :2.0; 

TF : 1 . 0 ; 

TFGUES: 1 .0; 

XGUES [1,1] : .3; 
XGUES [1,2] : 0 . 0 ; 
XGUES [2,1] : 0 . 0 ; 
XGUES [2,2] : —1 . 0 ; 
XGUES [3,1] : 2 . 0 ; 
XGUES [3, 2] : 2 . 0 ; 
NE : 8 ; 
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11.2. Sample output file 


The output (given on the following pages) of the state constraint example 
consists of the solutions for the states, costates, controls, and Hamiltonian for 2, 4, 
and 8 elements per phase. At the top of each page is the total elapsed computer time 
from the start of the program. On the two element case sheets is 10.40 secs. This 
is the time the code took to run the continuation method and the Newton-Raphson 
method for this case. This is a rather small number given the complexity of the 
problem and the fact that an accurate second-order Runge-Kutta method was used 
to solve the problem. The time at the top of the four element case is 11.94 which 
tells us that only 11.94 - 10.40 = 1.54 seconds was required to run the four-element 
case given the solution to the two element case. Finally, the desired eight element 
case solution was obtained in a total of 15.08 secs and only 3.14 secs from the four 
element case. Note that this time includes the extraction of nodal values and the 
production of the data files. This is a nonnegligible part of the total time. 

In summary, a third order state constraint problem which might have taken 
several days or weeks to program from scratch was solved in about 10 or 15 minutes 
with the general code. The simple input file is typed in a few minutes and a few 
minutes are required by MACSYMA to create the FORTRAN subroutines. After 
that, the program runs in a matter of seconds. 
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NODAL VALUES FOR THE STATES 


NUMBER OF ELEMENTS = 2 TOTAL ELAPSED TIME 


XI 

0 .OOOOOE+OO 
0.21250E+00 
0 . 30000E+00 
0.30000E+00 
0 . 21250E+00 
0.55511E-16 


X2 

0.10000E+01 
0 . 70000E+00 
0. 00000E+00 
0. 00000E+00 
- . 70000E+00 
- . 10000E+01 


X3 

0 . 20000E+01 
- . 44000E+01 
- . 12000E+01 
- . 12000E+01 
- . 44000E+01 
0 .20000E+01 


X4 

0. OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 


10.40 


TIME 

0 .OOOOOE+OO 
0 .25000E+00 
0 .50000E+00 
0 .50000E+00 
0 .75000E+00 
0 . 10000E+01 
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NODAL VALUES FOR THE STATES 


NUMBER 

OF ELEMENTS = 

4 TOTAL 

ELAPSED TIME = 

11.94 

XI 

X2 

X3 

X4 

TIME 

0 .OOOOOE+OO 
0 . 12385E+00 
0 .22506E+00 
0 . 28246E+00 
0 . 30000E+00 
0 . 30000E+00 
0 .28246E+00 
0 .22506E+00 
0 . 12385E+00 
- . 13878E-15 

0 .10000E+01 
0 . 98158E+00 
0 . 63780E+00 
0 .28062E+00 
0 . 277 56E-15 
0 . 27756E-15 
- .28062E+00 
- . 63780E+00 
- . 98158E+00 
- . 10000E+01 

0 .20000E+01 
- . 22947E+01 
- . 32057E+01 
- . 25091E+01 
- . 19809E+01 
- . 19809E+01 
- . 25091E+01 
- . 32057E+01 
- . 22947E+01 
0 .20000E+01 

0. OOOOOE+OO 
0. 00000E+00 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 

0 .OOOOOE+OO 
0 . 12500E+00 
0 .25000E+00 
0 .37500E+00 
0 . 50000E+00 
0 .50000E+00 
0 . 62500E+00 
0 .75000E+00 
0 . 87500E+00 
0 .10000E+01 
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NODAL VALUES FOR THE STATES 


NUMBER OF ELEMENTS - 8 TOTAL ELAPSED TIME - 


XI 

0 .00000E+00 
0 . 63949E-01 
0 . 12682E+00 
0 . 18238E+00 
0 .22738E+00 
0 . 2 608 6E+00 
0.28330E+00 
0.29595E+00 
0 . 30000E+00 
0 .30000E+00 
0 . 29595E+00 
0.28330E+00 
0 .26086E+00 
0 . 22738E+00 
0 . 18238E+00 
0.12682E+00 
0 » 63949E-01 
- . 15266E-15 


X2 

0 . 10000E+01 
0 . 104 64E+01 
0 . 96565E+00 
0 . 812 15E+00 
0 . 627 99E+00 
0 . 4 4308E+00 
0 . 27512E+00 
0 . 12965E+00 
- . 55511E-15 
- . 55511E-15 
- . 12965E+00 
- . 27512E+00 
- . 44308E+00 
- . 627 99E+00 
- . 81215E+00 
96565E+00 
- . 104 64E+01 
- . 10000E+01 


X3 

0 . 20000E+01 
- . 51643E+00 
- . 20665E+01 
- . 28453E+01 
- . 30479E+01 
- . 28695E+01 
- . 25051E+01 
- . 214 99E+01 
- . 19990E+01 
- . 19990E+01 
- . 21499E+01 
- . 25051E+01 
- . 28695E+01 
- . 30479E+01 
- . 28453E+01 
- . 20665E+01 
- . 51643E+00 
0 . 20000E+01 


X4 

0 . 00000E+00 
0 . 00000E+00 
0.00000E+00 
0 . 00000E+00 
0 . OOOOOE+OO 
0 . 00000E+00 
0. 00000E+00 
0. OOOOOE+OO 
0. 00000E+00 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 


15.08 


TIME 

0 .OOOOOE+OO 
0.62500E-01 
0 . 12500E+00 
0.18750E+00 
0 . 25000E+00 
0 . 31250E+00 
0 . 37500E+00 
0 . 43750E+00 
0 . 50000E+00 
0 . 50000E+00 
0 . 56250E+00 
0 . 62500E+00 
0.68750E+00 
0 . 7 5000E+00 
0 .81250E+00 
0 .87500E+00 
0 . 93750E+00 
0 . 10000E+01 
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NODAL VALUES FOR THE COSTATES 


NUMBER OF ELEMENTS - 2 TOTAL ELAPSED TIME = 


LI 

0.20480E+04 
0 . 20480E+04 
0 . 20480E+04 
- . 20480E+04 
- . 204 80E+04 
- . 204 80E+04 


L2 

0 . 66560E+03 
0 . 15360E+03 
- . 35840E+03 
- . 35840E+03 
0.15360E+03 
0 . 66560E+03 


L3 

0 .76800E+02 
- . 25600E+02 
- . 71054E-14 
- . 71054E-14 
0 . 25600E+02 
- . 76800E+02 


L4 

0 . OOOOOE+OO 
0 . 00000E+00 
0 .00000E+00 
0. 00000E+00 
0. 00000E+00 
0 .00000E+00 


10.40 


TIME 

0 .OOOOOE+OO 
0 . 25000E+00 
0 .50000E+00 
0 . 50000E+00 
0 .75000E+00 
0 .10000E+01 
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NODAL VALUES FOR THE COSTATES 


NUMBER OF ELEMENTS - 4 TOTAL ELAPSED TIME 


LI 

0 . 90935E+03 
0 . 90935E+03 
0 . 90935E+03 
0 . 90935E+03 
0 . 90935E+03 
- . 90935E+03 
- . 90935E+03 
- . 90935E+03 
- . 90935E+03 
- . 90935E+03 


L2 

0 . 33023E+03 
0 . 21656E+03 
0 . 10289E+03 
- . 10779E+02 
- . 124 45E+03 
- . 124 45E+03 
- . 10779E+02 
0 . 10289E+03 
0 . 21656E+03 
0 . 33023E+03 


L3 

0 . 51445E+02 
0 . 17271E+02 
- . 26947E+01 
- . 84517E+01 
0 . 46185E-13 
0 .46185E-13 
0 . 84517E+01 
0 . 26947E+01 
- . 17271E+02 
- . 51445E+02 


L4 

0 . 00000E+00 
0 . 00000E+00 
0 . 00000E+00 
0.00000E+00 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0. 00000E+00 
0. 00000E+00 
0. OOOOOE+OO 


11.94 


TIME 

0. 00000E+00 
0 . 12500E+00 
0 . 25000E+00 
0 .37500E+00 
0 .50000E+00 
0 .50000E+00 
0 . 62500E+00 
0 . 7 5000E+00 
0 .87500E+00 
0 . 10000E+01 


172 



NODAL VALUES FOR THE COSTATES 


NUMBER OF ELEMENTS = 8 TOTAL ELAPSED TIME 


LI 

0 . 79917E+03 
0.79917E+03 
0 .79917E+03 
0 . 79917E+03 
0.79917E+03 
0 . 79917E+03 
0 .79917E+03 
0 . 79917E+03 
0 . 79917E+03 
- . 79917E+03 
- . 79917E+03 
- . 79917E+03 
- . 79917E+03 
- . 79917E+03 
- . 79917E+03 
- . 79917E+03 
— . 79917E+03 
- . 79917E+03 


L2 

0 . 29734E+03 
0 . 24739E+03 
0 .19745E+03 
0 .14750E+03 
0 . 97549E+02 
0 .47601E+02 
- . 23476E+01 
- . 52296E+02 
- . 10224E+03 
- . 10224E+03 
- . 52296E+02 
- . 23476E+01 
0 .47601E+02 
0 . 97549E+02 
0 . 14750E+03 
0 .19745E+03 
0 .24739E+03 
0 . 2 9734E+03 


L3 

0 . 48774E+02 
0 .31751E+02 
0 . 17850E+02 
0 .70708E+01 
- . 58689E+00 
- . 51228E+01 
- . 65370E+01 
- . 48294E+01 
- . 26645E-14 
- . 26645E-14 
0 . 48294E+01 
0.65370E+01 
0 . 51228E+01 
0 . 58689E+00 
- . 70708E+01 
- . 17850E+02 
- . 31751E+02 
- . 48774E+02 


L4 

0 . 00000E+00 
0 .00000E+00 
0 . 00000E+00 
0 .OOOOOE+OO 
0 . OOOOOE+OO 
0 .00000E+00 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 


15.08 


TIME 

0. OOOOOE+OO 
0 . 62500E-01 
0 . 12500E+00 
0 .18750E+00 
0 .25000E+00 
0 . 31250E+00 
0 .37500E+00 
0.43750E+00 
0.50000E+00 
0 . 50000E+00 
0 . 56250E+00 
0 . 62500E+00 
0 . 68750E+00 
0 . 75000E+00 
0 .81250E+00 
0.87500E+00 
0.93750E+00 
0 .10000E+01 
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ALL VALUES FOR CONTROL AND HAMILTONIAN 


NUMBER OF ELEMENTS = 2 TOTAL ELAPSED TIME 


U1 

- . 76800E+02 
- . 25600E+02 
0.25600E+02 
0 . 128C0E+02 
0 . 15843E-13 
0.71054E-14 
- . 12800E+02 
- . 25600E+02 
0 . 25600E+02 
0 . 7 6800E+02 


U2 

0 . OOOOOE+OO 
0 . OOOOOE+OO 
0. 00000E+00 
0. 00000E+00 
0. 00000E+00 
0. OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 


U3 

0. OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 


HAMIL 

0 . 43008E+03 
0 . 92160E+03 
0 . 43008E+03 
0 . 92160E+03 
0 . 43008E+03 
0 . 43008E+03 
0 . 92160E+03 
0 . 43008E+03 
0 . 92160E+03 
0 . 43008E+03 


10.40 


TIME 

0 .OOOOOE+OO 
0 . 12500E+00 
0 .25000E+00 
0 .37500E+00 
0 . 50000E+00 
0 .50000E+00 
0 . 62500E+00 
0 . 75000E+00 
0 . 87500E+00 
0 . 10000E+01 
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ALL VALUES FOR CONTROL AND HAMILTONIAN 


NUMBER OF ELEMENTS - 4 TOTAL ELAPSED TIME 


U1 

- . 514 45E+02 
- . 34358E+02 
- . 17271E+02 
- . 72880E+01 
0 . 26947E+01 
0 . 55732E+01 
0 . 84517E+01 
0 . 42258E+01 
- . 58764E-13 
- . 46185E-13 
- . 42258E+01 
- . 84517E+01 
-.55732E+01 
- . 26947E+01 
0 . 72880E+01 
0 . 17271E+02 
0 . 34358E+02 
0 . 51445E+02 


U2 

0 . OOOOOE+OO 
0 .00000E+00 
0 .00000E+00 
0 .00000E+00 
0 .00000E+00 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. 00000E+00 
0 .00000E+00 
0. 00000E+00 
0. 00000E+00 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 


U3 

0. OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 


HAMIL 

0 . 24651E+03 
0 . 27045E+03 
0 . 24651E+03 
0 . 27045E+03 
0 . 24651E+03 
0 . 27045E+03 
0 . 24651E+03 
0 . 27045E+03 
0 . 24651E+03 
0 . 24651E+03 
0 . 27045E+03 
0 . 24651E+03 
0 . 27045E+03 
0 . 24651E+03 
0 . 27045E+03 
0 . 24651E+03 
0 . 27045E+03 
0 . 24651E+03 


11.94 


TIME 

0 .OOOOOE+OO 
0.62500E-01 
0 .12500E+00 
0.18750E+00 
0 . 25000E+00 
0 . 31250E+00 
0 . 37500E+00 
0 .43750E+00 
0 . 50000E+00 
0.50000E+00 
0 .56250E+00 
0 .62500E+00 
0 .68750E+00 
0 . 75000E+00 
0 . 81250E+00 
0 . 87500E+00 
0 . 93750E+00 
0 . 10000E+01 
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ALL VALUES FOR CONTROL AND HAMILTONIAN 


NUMBER OF ELEMENTS = 8 TOTAL ELAPSED TIME - 


U1 

-.4877 4E+02 
- . 4 02 63E+02 
- . 31751E+02 
- . 24801E+02 
- . 17850E+02 
- . 12461E+02 
- . 70708E+01 
- . 32419E+01 
0 . 58689E+00 
0.28549E+01 
0 . 51228E+01 
0 . 58299E+01 
0 . 65370E+01 
0 . 56832E+01 
0 .48294E+01 
0 . 24147E+01 
0 .44409E-15 
0 . 26645E-14 
- . 24147E+01 
- . 4 8294E+01 
- . 56832E+01 
- . 65370E+01 
- . 582 99E+01 
- . 51228E+01 
-.28549E+01 
- . 58689E+00 
0 . 324 19E+01 
0 . 70708E+01 
0 . 124 61E+02 
0 . 17850E+02 
0 .24801E+02 
0 .31751E+02 
0 . 40263E+02 
0 . 48774E+02 


U2 

0 .00000E+00 
0 .00000E+00 
0 .00000E+00 
0 .OOOOOE+OO 
0 . OOOOOE+OO 
0. 00000E+00 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0. 00000E+00 
0 .00000E+00 
0 .00000E+00 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0 . OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 


U3 

0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0 . OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0 . OOOOOE+OO 
0. OOOOOE+OO 
0 . OOOOOE+OO 


HAMIL 

0 . 20438E+03 
0 .20918E+03 
0 . 20438E+03 
0 . 20918E+03 
0 . 20438E+03 
0.20918E+03 
0.20438E+03 
0 . 20918E+03 
0 . 20438E+03 
0 . 20918E+03 
0.20438E+03 
0.20918E+03 
0 . 20438E+03 
0.20918E+03 
0.20438E+03 
0 . 20918E+03 
0.20438E+03 
0.20438E+03 
0 . 20918E+03 
0 . 20438E+03 
0 . 20918E+03 
0 . 20438E+03 
0 . 20918E+03 
0 . 20438E+03 
0 . 20918E+03 
0 . 20438E+03 
0.20918E+03 
0.20438E+03 
0 . 20918E+03 
0.20438E+03 
0.20918E+03 
0.20438E+03 
0.20918E+03 
0 .20438E+03 


15.08 


TIME 

0. OOOOOE+OO 
0.31250E-01 
0.62500E-01 
0.93750E-01 
0 . 12500E+00 
0 .15625E+00 
0 . 18750E+00 
0 .21875E+00 
0 . 25000E+00 
0 .28125E+00 
0 . 31250E+00 
0 . 34375E+00 
0.37500E+00 
0 . 40625E+00 
0 . 43750E+00 
0 . 46875E+00 
0 . 50000E+00 
0 . 50000E+00 
0 . 53125E+00 
0 .56250E+00 
0 .59375E+00 
0 .62500E+00 
0 . 65625E+00 
0 . 68750E+00 
0 . 71875E+00 
0 . 7 5000E+00 
0 . 78125E+00 
0 . 81250E+00 
0 .84375E+00 
0 .87500E+00 
0 . 90625E+00 
0 . 937 50E+00 
0 . 96875E+00 
0 .10000E+01 
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We conclude this chapter with a chart showing the actual clock time required 
to solve the above four problems. The four problems (to read the chart) are I, the 
free time trajectory optimization problem of Chapter 2, II, the control constraint 
problem of Chapter 6, III, the advanced launch vehicle of Chapter 8, and IV, the 
state constraint example problem presented above. 

The time to type in the different input files is approximately the same for all 
problems, around 5 to 10 minutes. The first column in Table 11.1 gives the clock 
time (in minutes) it took for MACSYMA to read in the input file and produce 
the two necessary FORTRAN subroutines to be used by the general code. These 
subroutines contain analytical first and second mixed partial derivatives for all the 
state equations. The second column gives the total elapsed time (in minutes) until 
the program generated all the requested results. This includes compilation and 
linking of the new subroutines created by MACSYMA. The third column gives the 
time required (in seconds) to generate each set of results for 2, 4, and 8 elements. 
These are the times given at the top of each output page as described above. It 
is seen by studying Table 11.1 that even the complicated advanced launch vehicle 
problem was solved in just 8 minutes and 20 seconds, from start to finish. Also, the 
majority of the time is used by MACSYMA to generate the subroutines. 
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Table 11.1: Setup and run times for the general code 


Problem 

Setup Time (minrsec) 

Total Time (min:sec) 

Run Time (sec) 2/4/8 

I 

5 : 05 

6 : 15 

7.80/9.28/11.14 

II 

2 : 45 

4 : 00 

3.68/4.80/6.12 

III 

6 : 10 

8 : 20 

21.76/25.60/33.38 

IV 

3:30 

4 : 40 

10.40/11.94/15.08 
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CHAPTER 12 


CONCLUSIONS AND FUTURE RESEARCH 


This thesis was concerned with a method for the solution of optimal control 
problems. The method, called the weak principle for optimal control, is based on 
time-domain finite elements. The goal of the weak principle was to develop an 
efficient and accurate algorithm in the hopes that it might be capable of produc- 
ing optimal trajectory solutions in a real-time environment. This would lead to 
tremendous savings in terms of time and money for many space related projects. 

Some of the characteristics and features of the weak principle are summarized 
below. 

(1) The necessary conditions of optimality are satisfied and all strong boundary 
conditions are transformed into weak boundary conditions. The weak princi- 
ple finds candidate extremal solutions, i.e., ones that satisfy all the necessary 
conditions. 

(2) Because all strong boundary conditions are cast in the form of natural or weak 
boundary conditions, then the same shape functions may be chosen for every 
optimal control problem. 

(3) The choice of shape functions allowed by the weak principle permits all inte- 
gration to be done by inspection, regardless of the degree of nonlinearity in the 
problem. Thus, no errors are introduced as the result of numerical quadrature. 

(4) One tremendous advantage of the integration being done by inspection is that 
algebraic equations may be derived prior to specifying the problem to be solved. 
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Because the form of the algebraic equations is known, it was possible to create 
a general algorithm for the solution of optimal control problems. 

(5) The algebraic equations which come from the weak principle possess a very 
sparse Jacobian. When this sparsity is exploited by way of a smart sparse 
matrix solver, a very efficient algorithm may be produced. 

(6) Symmetry in the solution is manifested by the weak principle. Problems, such 
as those presented in Chapter 2, that possess a state, costate, or control that 
is analytically symmetric, will also have symmetric approximated solutions 
obtained by the weak principle. 

These features of the weak principle make it a powerful and versatile tool for 
solving optimal control problems. Particular attention was given to the applica- 
tion of the method to advanced launch vehicle guidance. A real-time algorithm 
was to be developed which requires that the optimal trajectory be computed in a 
small time interval as compared to the total time interval. This thesis presented 
a solution to this optimization problem which was obtained in about 5.5 secs as 
compared to the 360 second time span of the entire trajectory. This fact, coupled 
with the accuracy demonstrated by the weak principle, make real-time trajectory 
optimization a possibility. 

A general code for the solution of optimal control problems was developed 
based on the weak principle. It was possible to create a general purpose, robust, and 
efficient algorithm because the algebraic equations are known before the problem is 
specified. The most promising feature of the general code is the reduction in setup 
time for any given problem over any other existing formulation. Problems that 
might take weeks to program and solve can be solved in a matter of minutes with 
the general code. This can lead to tremendous savings in terms of time, money, and 
human resources. 

There is still work that can be done for developing the weak principle further. 
Below are a few ideas. 
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(1) The algebraic equations are derived in terms of the midpoint values of the shape 
functions. For problems with staging or state constraints, certain nodal values 
appear and must be recovered around the unknown time(s). It may be easier 
if all the algebraic equations were rewritten in terms of the nodal values. This 
may also be more efficient since the nodal values need to be recovered anyway. 
However, the resulting Jacobian may not be as sparse as with the midpoint 
values. 

(2) It was noted that the algebraic equations were all linear in terms of the costates. 
It may be possible to use MACSYMA to symbolically solve for the costates in 
terms of the other unknowns. The advantages of this would be two-fold. One 
is that the number of equations would be cut almost in half and there would 
be no associated decrease in the percentage of zeros in the Jacobian. Second, 
there would be no need to obtain initial guesses for the costates. This is helpful 
because there is often times little or no physical insight into the magnitude of 
the costate variables. 

(3) It may be possible to generate a canonical form for constraints that will allow 
for easier programming of multiple constraints. A canonical form, which is 
successfully used in [39], could have several advantages, including that it might 
allow for multiple constraints and that there may also be a way of eliminating 
the users need for estimating a priori and iterating on the number of phases 
present in the solution. 

(4) Work can be continued on the weak principle to remove the singularity asso- 
ciated with state constraint problems when using the necessary conditions of 
Ref. 65. This would also allow more uniformity in the programming. 

(5) The weak principle could be extended to include singular control problems 
(see [1]). Singular control problems are characterized by the control appearing 
linearly in the formulation. Since the optimality condition does not determine 
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the control, successive time derivatives are taken until the control does appear. 
It should be possible to incorporate this type of problem into the weak principle. 

(6) Finally, the general code can be continually updated. Along with the sugges- 
tions above that can be incorporated into the general code, there axe other 
potential savings to be explored. One is the generation of the subroutines by 
MACSYMA. Perhaps there are more efficient ways of producing the needed 
code, or perhaps Mathematica (another symbolic manipulator) would be bet- 
ter. In addition, the general code would be very useful if it were available for 
desktop computers. One drawback of existing codes is that many of them only 
r un on large, and sometimes inaccessible, machines. A fast and efficient code 
to run on desktop computers would be of great value to industry and academia. 
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APPENDIX A 


HAMILTON’S WEAK PRINCIPLE FOR DYNAMICS 

The potential of obtaining a direct solution in the time domain is very much 
analogous to obtaining the solution of a beam deflection problem with the beam 
axial coordinate broken into several segments or finite elements. In the present case, 
however, it is the time interval which is broken into segments; thus, the phrase “finite 
elements in time” has been adopted by several investigators. 

Only recently has a mixed formulation of Hamilton’s Weak Principle (HWP) 
been investigated as a computational tool for finite elements in time [47]. In this 
section, the mixed form of HWP is derived and its application to dynamics problems 
is illustrated. 

A.l. General Development 

To this aim, let us consider an arbitrary holonomic mechanical system. The 
configuration is completely defined by a set of generalized coordinates q. Further, 
let us denote with L(q , q, t) the Lagrangean of the system, Q the set of nonconserva- 
tive generalized forces applied to the system, and p = dLj dq the set of generalized 
momenta. The generalized coordinates q should be piecewise differentiable and the 
generalized momenta p will have discrete values at to and tf. (For a more math- 
ematically rigorous discussion, see [53].) Then the following variational equation, 
known as HWP [45], describes the real motion of the system between the two known 
times to and tf. 
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(A.l-1) 


f 6Ldt+ f 8q T Q dt = Sqjfif — Sq^po 

where 8q, the variation of q , should be of the same class of functions as is q, 8qf = 
8q(t = tj ) and <$90 = 8q(t = to)- Although HWP contains p in the form of discrete 
values at the end points denoted by the hatted quantities, this particular variational 
equation is said to be in displacement form because it only involves the variation 
of q. Keeping p distinct from dL/dq allows for accurate extraction of the momenta 
without differentiation of q. Although this formulation has been shown to be of 
practical use in dynamics (see [45] and [46]), an even more useful formulation may be 
derived if independent variations in both displacements and momenta are allowed, 
resulting in a mixed formulation. 

In order to derive the mixed formulation, let the Hamiltonian be defined as 

H{q,p, t) = p T q - L(q, q, t ) (A.l - 2) 

Taking the variation of Eq. (A. 1-2) and substituting for 8L in Eq. (A.l-1) results 
in 



q + 8q T p — 8H + 8q T Q ) dt = 8qjpf — dqjpo 


Now, introducing 


(A.l -3) 


and 




= lim 
<— tt 


5f(0 


and q\tf — lhn ?(<) 


?k=9(* 0 ) and q\t,=q(tf ) 


( 12 - 1 ) 


( 12 - 2 ) 


184 

C'3 



then continuity between the values of q and p on the interior and q and p on the 
boundary is weakly enforced by adjoining 8a T (q - q)\ l / 0 to Eq. (A. 1-3) where 8a 
is a set of discrete unknown Lagrange multipliers defined only at to and tf. The 
resulting equation is 


f {8p T q + 8q T p — 8H + 8q T Q) dt = 8qjpf — 8q£p 0 + 8a T (q g)| t g (A.l 4) 

Jt o 

It is now possible to choose 8a = 8\ without changing any necessary conditions. 
Also, to finish the development, the first term in Eq. (A. 1-4) is integrated by parts 
yielding 

/ (8q T p- 8p T q - 8H + 8q T Q) dt = 8q 1 fpf - 8qoPo - 8p T jqf + 8p^qo (A.l - 5) 
Jt o 

This is called a mixed formulation because it contains independent variations of q 
and p . It is also in the ^^eakest” possible form in the sense that all boundary condi- 
tions are of the natural type, enforced by the variational equation for unconstrained 
variations. Note that now p aud q should have discrete values at to and tf^p and 
q should be piecewise continuous, and 8p and 8q should be piecewise differentiable 

(O' 1 ). 

There are two main advantages of the mixed formulation over the displacement 
formulation. The first advantage is that the mixed formulation generally provides 
a more accurate solution for a given level of computational effort than does the 
displacement formulation. The second advantage is that a simpler choice of shape 
functions is allowed. Note in Eq. (A.l— 5) that time derivatives of 8q and 8p are 
present. However, no time derivatives of q and p exist. Therefore, it is possible to 
implement linear shape functions for 8q and 8p and constant shape functions for q 
and p within each element. 
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A. 2. Finite Elements in Time 


Let the time interval from to to tf be broken into N equally spaced elements. 
The nodal values of these elements are <, for i = 1,. . . ,N + 1 where to = t\ and 
tf = t N+i- A nondimensional elemental time r is defined as 


_ t - U _ t-t, 
t*-f" i tj At * 


(A. 2 — 1) 


The linear shape functions for 


the virtual coordinates and momenta are 


8q = Sq,( 1 - r) + Sq i+1 T 

Sp = 6pi( 1 -T) + Sp i+ iT 


(A. 2 — 2) 


For the generalized coordinates and momenta 


and 


<1 = 


p = 



if r = 0; 

qi 

if 0 < r < 

qi+ 1 

if r = 1 

Pi 

II 

O 

Pi 

if 0 < t < 

Pi+i 

if t = 1 


(A. 2 — 3) 


(A. 2 — 4) 


It is important to understand that qi = g(<o),Pi = p(to),qN+i = and pn+i - 

p(t/). In other words, the hatted values of q and p at the beginning and end of our 
time marching scheme are the discrete values of q and p that are needed in the mixed 
formulation. When these shape functions are substituted into Eq. (A. 1-5), one can 
either generate an implicit time-marching procedure for nonlinear problems or apply 
standard finite element assembly procedures to solve periodic or two-point boundary 
value problems [48]. When this formulation is applied to the linear oscillator, a 


186 



time-marching algorithm emerges that is unconditionally stable [45]. Higher-order 
(so-called p- version) elements could be developed [69], and they would certainly be 
attractive for linear problems or for nonlinear problems with nonlinearities of low 
order. For nonlinear problems in general, use of the crude shape functions allowable 
with the mixed method would seem to be more efficient than use of higher-order 
shape functions in a p- version. The reason for this is that, with the exception of the 
term involving Q, which may contain time explicitly, all element quadrature can be 
done by inspection regardless of the order of the nonlinearities. 

A. 3. Example: A Nonlinear Initial- Value Problem 

Applying the shape functions of Eqs. (A. 2-2 - A.2-4) to Eq. (A. 1-5) for an 
initial value problem, a recursive set of nonlinear algebraic equations is obtained of 
the form 


fj(qi,qi+uPi,Pi+ 1) = 0 i = (A. 3 — 1) 

where n is four times the number of degrees of freedom of the system. Eq. (A. 3-1) 
can be solved by a Newton- Raphson method yielding an implicit time-marching 
procedure. The key advantage of using finite elements and a weak variational ap- 
proach over numerical integration is that the solution (for linear problems) is stable 
for all time steps. In other words, no matter how large a time step is used, a finite 
approximation of the solution will be obtained. This unconditional stability is ob- 
tained without ad hoc procedures such as selective or reduced element quadrature 
which are necessary in displacement formulations. 

Also noteworthy of the finite element discretization is that the midpoint values 
of q and p are just the average values of the adjoining nodal values, or 


2 qi =qi + qi + 1 

2 pi =Pi + Pi+i 


(A. 3 — 2) 
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Thus, it is possible to cut the number of equations and unknowns in half. This can 
be very useful for a multi-degree of freedom problem in terms of computer savings. 

Consider a simple pendulum composed of a lumped mass m and a weightless 
bar of length £ (see Fig. A.l). The single generalized coordinate q is the angular 
displacement of the bar from the vertical. Denoting the kinetic energy of the system 
with K and the potential energy with V, then we may define the following: 


K = \mPq 2 

2 y 

V = mg£( 1 — cos <7) 

L — I\ -V 

dL 02- 
p = — = mi q 

dq 

P* 

H = + mg£(l-cosq) 


(A. 3 — 3) 


There are no nonconservative forces Q applied to this system. 

Substituting t = t { + rAt; from Eq. (A.2-1), along with Eq. (A.3-3), and 
substituting the shape functions defined in Eqs. (A.2-2 - A. 2-4) into Eq. (A. 1-5) 
we obtain (for i = 1,2 ,N) 


A t. 



A U ) 


f tpi+i - hi \ 

V At, ) 


Pi - mgt sin qi [%(1 - r) + Sq i+1 r] 

~ - r ) + hi+ir] j dr 


~ ^Qi+iPi+l + Sp i+1 q i+1 + Sqipi - Spiqi = 0 


(A.3 - 4) 


Carrying out the integration by inspection and setting the coefficient of each virtual 
quantity (Sqi, Spi , £5,4-1 , and £p,4-i ) to zero, the following four independent equations 
for each value of i are obtained. 
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* _ mgPAtisinqi 

P'~P' 2 

mgPAti sin qi 
Pi ~Pi + 1 2 

PiAti 
qi qi 2 mi 2 

piAti 

q ' +1 qi 2 mP 


= 0 
= 0 
= 0 
= 0 


(A. 3 — 5) 


There are six unknowns; however, for an initial-value problem, we will specify qi 
and pi and solve for the remaining unknowns as outlined below. Thus, Eq. (A. 3-5) 
is of the form of Eq. (A. 3-1). 

Recall that i ranges from 1 to N. To start with, i — 1 and <71 and pi (i.e. the 
initial conditions) are specified. Now, solve for <7 i,Pi>< 72, and p 2 - Next, let i = 2, 
use the known <72 and f >2 (we just found those values), and solve for the new four 
unknowns. This process is repeated until i = N. 

For this simple pendulum example, the variables will be nondimensionalized as 
follows. If we define u> 2 = g/£, then a dimensionless time step At may be defined 
that does not vary with i so that At = u>Ati. Also, instead of solving directly for 
p, the dimensionless p/mPu) will be solved for. 

The initial conditions of the pendulum are <71 = 60° and pi = 0.0. The equa- 
tions will be solved for At = 0.4, 0.8, and 1.6. Graphs of the solutions are shown 
in Fig. A. 2 and Fig. A. 3 and compared to the exact elliptic integral solution [78]. 
(Since the element midpoint values are simply the average of the element nodal 
values, only the nodal values are given for these results.) From Figs. A. 2 and A. 3, 
it is easily seen that At = 0.4 gives acceptable results for both displacement and 
angular momentum. Also, note that even the large 1.6 time step yields a finite 
approximation of the exact solution. 
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Fig. A.l: Nomenclature for example 

A simple pendulum composed of a lumped mass m and a weightless bar of length l 
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Dimensionless Momentum 



Dimensionless Time 


Fig. A. 3: Dimensionless momentum p/m£ 2 u versus dimensionless time 
Results for three values of the time step At and the exact elliptic integral solution 
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APPENDIX B 


INITIAL-VALUE ODE’s 

A very useful idea came from the work described in Chapter 2 and Appendix A 
concerning the solution of first-order ordinary differential equations. Of course, 
initial-value solvers are ubiquitous, so the goal of this work was to produce an 
easier and more efficient method to solve the differential equations. Furthermore, 
it was desired to employ a symbolic environment (such as MACSYMA) to perform 
all the necessary computations and write certain FORTRAN subroutines. This 
is similar to what is done by the general code described in Chapter 11, only this 
program is much simpler. This aspect of the work fully automates the procedure of 
solving initial-value problems and minimizes the amount of user interaction. 

B.l. General Development 

The problem may be stated as follows: Given a set of n first-order ordinary 
differential equations of the form 

x = f(x,t ) t\ < t < t 2 (B.l — 1) 

with x(to) specified, find x(t). 

To begin, a weighted residual method is used and Eq. (B.l-1) is integrated 
over the time interval of interest. 



Now, as was done in Chapters 2 and Appendix A, the strong boundary condi- 
tions are transformed to natural boundary conditions by weakly enforcing continu- 
ity of the states at the initial and final times. This is done via a discrete Lagrange 
multiplier which we identify as 6 A. Eq. (B.l-2) now becomes 

f 6X t [x — dt + SX T (x — x)|** = 0 

J <i 

Integrating the above equation by parts results in 

(a T x + £A r /) dt - SX T x\ t t l= 0 (B.l - 4) 

Once again introducing a nondimensional time r as 

(B.l -5) 

a linear shape function for SX of the form 


t-t i _t-ti 

T 2 - <i “ At 



(B.l -3) 


<5A = <$Ai(l — r) + ^A,-|-it 

and a piecewise-constant shape function for x of the form 


(B.l -6) 


Xi 

1—* • 
II 

i? 

x t 

if 0 < r < 1; 

^.+1 

if r = 1 


(B.l -7) 


then these values are substituted into Eq. (B.l-4). Carrying out the integration, 
the following set of algebraic equations is obtained 


- At f ' 

Xi - —fi = 

At , . 

Xi H" 2 /« — ®»+i 


(B.l -8) 


194 



where /, = f(x = Xi,t = ti). 

Eq. (B.l-8) yields a time-marching algorithm since the value of A ti is specified 
at each time step. In practice, one solves the first of Eq. (B.l-8) for Xi, and then 
obtains the nodal value x,+i from 

= 2 Xi — Xi (B.l — 9) 

This process then repeats until the final desired value of x is reached. 

B.2. MACSYMA: A Symbolic Manipulator 

Eq. (B.l-8) is a set of very simple nonlinear equations. These equations can 
be solved by a Newton- Raphson method. Since we are time-marching, the previous 
known nodal values serve as initial guesses. To date, the roots of the equations have 
always been found using these initial guesses. 

Using the above method, it was necessary to write a new code (or a large part 
of it) incorporating the new equations for each initial-value problem. Also, it was 
necessary to take explicit derivatives of the /’ s and code those values in for the 
linearization process. Therefore, it was desirable to combine FORTRAN code and 
MACSYMA code to automate our initial- value solver. 

MACSYMA is a large symbolic manipulation program developed at the MIT 
Laboratory for Computer Science [77]. This program has many capabilities which 
include taking an analytical expression, finding an analytical derivative, writing 
FORTRAN code for the analytical expression, and splicing the FORTRAN code 
into a subroutine template file. This is exactly what we needed MACSYMA to do. 

The general procedure to solve initial- value ordinary differential equations con- 
sists of three batch files, two subroutines, and a main program. Each of the above 
six components of the procedure is briefly described below. The whole procedure is 
very compact and runs very efficiently. 
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One batch file, “run. bat,” is the top-level command file that controls the whole 
procedure. Its commands include putting the user into an input file in order to 
input the /’ s for MACSYMA to read. Next, run.bat staxts up MACSYMA. At this 
point other batch files take over. Upon return to run.bat, the commands include 
compiling the MACSYMA written subroutine, adding the subroutine to an archive 
library, compiling the main program with the updated archive library, and finally 
running the program. 

The MACSYMA input file, named “input. macsyma,” prompts the user to enter 
the order of the system, n, and then the /’ s of the system in MACSYMA form. 
For example, if one wishes to solve the scalar system x = x, then “input. macsyma” 
would contain the two lines “N:l;” and “G[l]:x(l);’\ 

Once MACSYMA is started by “run.bat,” the user must type “batch(run);”. 
This batch file calls another batch file which loads the user supplied equations. Then 
MACSYMA is asked to evaluate the expressions in the template file and splice in 
the FORTRAN code. Finally, MACSYMA is terminated and control is returned to 
“run.bat.” 

The template file is very short and simple. It is nothing more than a FOR- 
TRAN subroutine (IVMAC) with MACSYMA expressions put in some places. After 
MACSYMA is called, all MACSYMA expressions in the template file are replaced 
with legal FORTRAN statements. This interaction between FORTRAN and MAC- 
SYMA was the key to automating the procedure and minimizing user interaction. 

The largest file, containing only 83 lines of code, is the subroutine SOLVER. 
This routine linearizes the discretized algebraic equations, calls IVMAC for values 
of / and derivatives of /, calls some Harwell subroutines to solve the equations, and 
writes the data to an output file. 

Finally, the main program, IVODE, asks the user for initial conditions, the 
time interval, and the time step to take. IVODE calls SOLVER and computes the 
elapsed computer time. 
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The above work has been a key part of the research effort. Although the initial- 
value problem is much easier than the optimal-control problem, the work indicated 
that a general procedure to solve optimal-control problems is realizable. The above 
outlined procedure must simply be broadened to handle the more complicated equa- 
tions that come from optimal-control problems. 
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APPENDIX C 


APPLICATION TO BEAM THEORY 

The weak principle for optimal control developed in this thesis and, in par- 
ticular, the general code described in Chapter 11 may be used to solve virtually 
any problem that can be cast in the proper form. Specifically, if one can identify 
a performance index and state equations, then the general code could be used to 
solve these problems. Problems from areas such as chemical engineering, robotics, 
electrical engineering and elasticity could be solved. This appendix examines how 
simple beam problems can be solved with the general code of Chapter 11. 

C.l. Transformation of Beam Problem 


Consider a simple cantilever beam of length L with a distributed load q, an 
end load Fl, and an end moment Ml- The deflection and slope of the neutral axis 
of the beam are signified by v and ft respectively; and the curvature is denoted by 

K. 

The first thing to do is identify elements of the beam problem with elements 
of the optimal control problem. From simple elasticity, we know that the first 
derivative of deflection with respect to the beam axis yields the slope, and the 
second derivative gives the curvature. These two equations can be written as two 
first-order equations as 


v' — ft and ft' = k 


(C.l - 1) 
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It is easy to identify v and ,8 as the two states and k as the control variable. What 
remains is to identify a suitable performance index. It is well known in elasticity 
that the variation of the strain energy equals the variation of the work produced by 
external forces. Therefore, a performance index may be written as 

( — qv) dx - Flvl ~ MlPl (C.l — 2) 

where El is the bending stiffness of the beam and x measures the distance along 
the axis of the beam. The Hamiltonian can now be identified as 



H = 


EIk 2 

2 


— qv + A v f3 + 


(C.l — 3) 


The costate equations yield additional insight into the problem. The equations are 


\' v = -— = q and A'^ = -— = -A„ (C.l -4) 

From these differential equations, we identify A„ = — F and A^ = — M, or in 
words, the costates yield the shear force and moment distribution along the beam. 

The beam problem has now been transformed to the equivalent optimal control 
problem and is ready for solution. 

C.2. Example 

A simple example problem is now considered. Consider a cantilevered beam 
with a tip load. Let El = 10 6 psi, L = 100 in, and the tip load Fi = 30 lb. The 
states will have the following constraints at the final time 


ipi=v- 


F l L 3 

3EI 


and 


4>2 = 0 - 


f l l 2 

3EI 


(C.2 - 1) 
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These constraints impose that the beam have a deflection equal to that which it 
would have with no constraint, however, the slope of the beam at the end is changed 
so that is lies along a line joining the origin and the point of deflection at the end 
of the beam. 

The equations and boundary conditions were put into the general code (as 
described in Chapter 11) and the solution was readily found. Figs. C.l - C.4 show 
the deflection, slope, shear force, and moment distributions for 2 and 4 elements 
and the exact answer. An 8 element case is also given for the shear force in Fig. C.3. 
In all the graphs, the accuracy of the finite element method is once again seen. 
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Fig. C.l: Deflection versus axial coordinate 
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Slope (rad) 



Fig. C.2: Slope versus axial coordinate 
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Shear Force 
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Fig. C.3: Shear force versus axial coordinate 
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